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PREFACE 


This report (Volume I of three volumes, with Volume II to be NASA CR-141743 
and Volume III to be NASA CR-141744) partially documents the development of a com- 
puter program (NATA) for calculating the flow in arc-heated wind tunnels and the 
conditions on models tested in such reentry simulation facilities. The objective 
was to provide a means for predicting and interpreting test conditions used in 
the experimental evaluation of thermal protection materials for reentry vehicles 
such as the Space Shuttle Orbiter. Much effort was expended to make the program 
reliable and easily usable by engineers without great expertise in gas phase 
chemical kinetics, gas transport properties, thermochemistry, and computer pro- 
gramming. The capabilities of NATA are summarized concisely in the abstract, 
and in somewhat more detail in the Introduction (Section 1) . Comparisons with 
available experimental data indicate that the results produced by the program 
have a useful level of accuracy. 
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MATHEMATICAL SYMBOLS 


All mathematical symbols used in this report are defined in the text where 
they are first used. In addition, the definitions of symbols used in several 
different sections of the report are summarized below for convenience of 
reference. 


Latin Symbols 


A g* 


Effective area ratio, A' e /A' e , 

Effective cross sectional area of the inviscid flow at a given axial 
position 

Effective cross sectional area of the inviscid flow at the sonic point 

Effective area ratio in the non-equilibrium solution by the inverse 
method 


A g 

A g0 

A: 


A ki 


"P 

‘ri 

S 

C pe 

Si 

D am 

E 

E 0 ; 


Geometric area ratio , A ' / A ' 0 

Geometric cross sectional area of a nozzle or channel at a given axial 
position 

Geometric cross sectional area at the throat 

Constant coefficient in formula for forward reaction rate constant 
for the i th reaction 

Inverse matrix of the matrix a~ 

Number of chemical elements present in the gas mixture 
Specific heat of the gas at constant pressure 
Specific heat for the j f b species at constant pressure 
Molar heat capacity of the gas at constant pressure 
Molar heat capacity for the electrons 

Molar heat capacity for the j ^ 1 species at constant pressure 
Atom-molecule binary diffusion coefficient 
Molar internal energy 

Activation energy in formula for forward reaction rate constant for 
the 1 th reaction 


Energy of electronic excitation for the state in a molecule 
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K e Equilibrium thermal conductivity 

Kj Chemically frozen thermal conductivity 

K. Equilibrium constant for the i 1 * 1 reaction 

m Mass flux, pa 

m e Electron mass 
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N pr Prandtl number, c p p/K 
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Pressure 

Partial pressure of the j th species 
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Rate factor in equation for rate of change of the concentration of a 
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species due to the i cn reaction, (p /u)kfj n yjj 

Rate of energy addition to the electron gas per unit volume 

Mean radiative energy loss from the gas inN Q occurrences of the 
reaction in the forward direction 

Radiative power loss from the gas per unit volume 

Negative of the mean radiative energy loss from the gas in N 0 occurrences 
of the i* reaction in the reverse direction 

Stagnation point heat flux 

Heat flux to the wall 

Number of gram-atoms of the k th element per mole of the cold gas mixture 

Number of reactions in the reaction system 
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MATHEMATICAL SYMBOLS (Cont'd) 


T e Absolute electron temperature 
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u Flow velocity 

V Molar volume 
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THE NATA CODE - THEORY AND ANALYSIS 


By W. L. Bade and J. M. Yos 
Avco Systems Division 
Wilmington, Massachusetts 


1. INTRODUCTION 


The NATA* code was developed by Avco Systems Division from 1968 through 
197.4, under the sponsorship of the NASA Johnson Space Center, to meet a need for 
predicting and interpreting the flow conditions in electric-arc heated wind tun- 
nels. Test facilities of this type are used at the Johnson Space Center, and 
elsewhere, for evaluating thermal protection materials for reentry systems. The 
present report documents the theory and analysis upon which the current version 
of NATA is based. A user's manual and a programmer's manual for the code are to 
be issued separately. 

The earliest version of NATA was based upon a computer program developed at 
Cornell Aeronautical Laboratory (CAL) (ref. 1). This CAL program performed cal- 
culations of quasi-one-dimensional flow through a nozzle of specified geometry. 
The flow was assumed to start from a state of thermochemical equilibrium at spe- 
cified temperature and pressure in an upstream reservoir. The program included 
options for flow solutions based on chemical equilibrium, frozen chemistry, and 
chemical non-equilibrium with specified reaction rate constants. Boundary layer 
effects were neglected. There was no provision for calculating conditions on 
models immersed in the flow. The program used fixed-format input; i.e., the 
numerical data were punched, without alphanumeric identifiers, in prescribed 
fields of the input cards. The input was voluminous (100 to 200 cards per case), 
as all of the gas species properties and reaction rate data had to be read in. 

The output was in the form of non-dimensional ratios of the gas flow properties 
to standard values. 

During the development of NATA, the basic calculations of frozen, equilib- 
rium, and non-equilibrium flow as coded in the CAL program have been retained 
with only minor changes. Many new features have been added, including a laminar 
boundary layer calculation, internal gas transport property calculations, options 
for specifying the reservoir conditions by input of the total mass flow together 
with the reservoir pressure or the stagnation enthalpy, and computations of the 
heating and stresses on models immersed in the flow. In addition, a radical 
revision of the input and output arrangements has been carried out to make the 
code easier to use. The fixed-format input of the CAL program has been replaced 
by namelist input, in which only the variables whose values are to be changed 
have to be read in. The control inputs are all preset to standard values which 
are generally satisfactory. All species and reaction data for several standard 
gas models, including argon-free air, are compiled into the program and can be 
invoked by specifying a single input variable. Curvefits to the geometries of 
several of the nozzles and channels used at NASA Johnson Space Center are also 
included in the program. As a result of these changes, the number of input cards 
required per case has been reduced to 4 to 8 (typically) . The output has been 


•Acronym for Non-equilibrium Arc Tunnel Analysis. 
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expanded to include additional quantities, such as the boundary layer properties, 
and is now presented in dimensional form to allow direct comparison of the code 
results with experimental data. Each output variable is labeled with an alpha- 
numeric identifier. In addition, much effort has been expended to ensure that 
NATA will almost always generate a successful solution when run with the standard 
values of the control inputs. A major objective of these revisions was to make 
the code usable by gas dynamicists who are not necessarily expert in thermo- 
chemistry, transport property theory, chemical kinetics, or computer programming. 

The remaining parts of this introduction outline the capabilities of NATA, 
discuss its limitations, and describe the subsequent sections of the present 
report. 


1.1 Programming 

NATA is a Fortran IV program consisting of a main program and 68 subroutines. 
The source deck contains approximately 8500 cards . The program exists in two 
versions, one for use on the IBM 360 system, the other for the Univac 1108. The 
IBM 360 version is entirely in double precision, whereas the Univac version is a 
single precision program with some double precision arrays and subroutines. One 
version can be converted into the other by inserting or removing the cards which 
type all floating-point variables as double precision in the IBM 360 version. 

The IBM 360 version requires about 410K bytes of core storage, including 
buffers. The Univac version is run on the 1108 using overlay, and fits into the 
two-bank processors at NASA/ JSC with about 2000 words of storage to spare. 


1.2 Gas Models 

NATA provides four options for specifying the species properties and re- 
actions defining the gas model. 


1.2.1 Standard Gas Models 

Six compiled- in standard gas models are available in the code. One of these 
is a model for argon-free air at low and moderate temperatures where N0+ is the 
only important ion species. Another is a high-temperature air model including 
five atomic and molecular ions. There are two planetary atmosphere models for 
mixtures of 75 mole percent CO 2 , 20 mole percent Ar, and 5 mole percent N 2 . 
Finally, there is an electronic non-equilibrium (two-temperature) model for argon 
and one for helium. These rare-gas models include a non-equilibrium treatment 
of electronic excited states. Any of these standard models can be selected by 
input of a single index value. 

1.2.2 Standard Gas Models with Input Mole Fractions 

The elemental composition of the gas mixture is specified in terms of the 
mole fractions of "cold species" which are stable at low temperature. For exam- 
ple, the weight fractions of C, 0, N, and Ar in the planetary atmosphere models 
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are determined from the assumed mole fractions of CO 2 , N2> and Ar in the gas be- 
fore it is heated. Optionally, these cold species mole fractions can be set in 
the input without disturbing the standard species properties or the reaction sys- 
tem. Thus, the standard models for air and the planetary atmosphere can easily 
be applied to mixtures of the same cold species in different proportions. 


1.2.3 User-Generated Gas Models with Standard 
Species and Reactions 

A gas model based upon the available, compiled-in species and reactions, 
but with a non-standard selection of species and reactions, can be set up in the 
input. For example, some of the reactions can be deleted from a standard model 
to assess their effect upon the solution of a particular case. Input of such a 
non-standard gas model requires about 4 to 8 input cards, depending upon the 
numbers of species and reactions included. 


1.2.4 Input of Species and Reaction Data 

Finally, the basic data for species and reactions can be set in the input. 
The data for standard species and reactions can be changed in part or in whole, 
and new species and reactions can be introduced. About two input cards are re- 
quired per new reaction, and about 3 to 4 cards per species. 


1.3 Transport Properties 

NATA includes a self-contained capability for computing the transport 
properties of gas mixtures, based upon the temperature, pressure, and mole frac- 
tions given by the equilibrium or non-equilibrium calculations of the flow solu- 
tion. The transport collision cross sections for the standard species are 
provided by compiled-in data. Collision cross section data can also be specified 
in the input. The transport properties calculated are the viscosity, the chemi- 
cally frozen Prandtl number, the atom-molecule Lewis number, and the electrical 
conductivity. The viscosity and Prandtl number are used in the calculations of 
the boundary layer on the nozzle wall and in calculations of stagnation-point 
heating on models and of the heat transfer to wedge models. The Lewis number is 
used in stagnation-point heat transfer calculations. The electrical conductivity 
is not used in NATA, but is printed out. 

1.4 Flow Geometry 

The geometries of nozzles and channels are specified, in NATA, by means of 
curvefits to their profiles. A profile is the curve of intersection of the 
inner surface of a nozzle or channel with a symmetry plane. For an axisymmetric 
nozzle, there is one such profile. For a rectangular channel, there are two 
profiles. Each profile is represented as a sequence of straight lines and cir- 
cular arcs, joined end-to-end with value and slope continuity. A separate 
computer program (NOZFIT) is available for generating such curvefits from 
nozzle design data. 
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If boundary layer displacement effects are neglected, the inviscid flow de- 
pends only upon the ratio of the nozzle cross sectional area to the area at the 
throat, not upon the shape of the cross section. However, the convergence or 
divergence of streamlines in the boundary layer does depend upon the type of 
nozzle geometry, and affects the rate of boundary layer growth. To take such 
boundary layer effects into account, NATA contains explicit treatments of three 
types of nozzle geometry: axisymmetric nozzles, two-dimensional nozzles, and 

rectangular channels. 

1.5 Reservoir Conditions 

The flow is always assumed to start from a state of thermochemical equilib- 
rium in an upstream reservoir. For a gas of given elemental composition, the 
thermochemical state is a function of two variables, and its specification there- 
fore requires two inputs. NATA provides three options for the input specification 
of reservoir conditions. 


1.5.1 Temperature and Pressure 

The basic method is the one used in the original CAL program (ref. 1), 
namely, direct input of the reservoir temperature and reservoir pressure. NATA 
contains a subroutine for determining the equilibrium species mole fractions 
from these data and the species thermochemical properties. 


1.5.2 Pressure and Mass Flow 

The upstream stagnation pressure can be measured easily in arc heated wind 
tunnels in which the flow stagnates in a plenum chamber downstream of the arc 
heater and upstream of the throat. However, the reservoir temperature is diffi- 
cult to determine. For this reason, an option is provided to calculate the 
reservoir conditions from data on the reservoir pressure and the total mass flow 
(which is easily measured) . When this option is used, the reservoir temperature 
is estimated and the resulting total mass flow is calculated based on an equilib- 
rium flow solution from the reservoir to the throat. An iteration is then 
carried out to determine the reservoir temperature corresponding to the input 
mass flow. If the boundary layer is to be included in the main flow solution, a 
correction is made for the effect of the displacement thickness at the throat. 


1.5.3 Mass Flow and Stagnation Enthalpy 

Many arc heated wind tunnels lack an upstream stagnation region. In such 
facilities upstream pressure measurements do not give the effective reservoir 
pressure. To deal with such cases, NATA provides a third option in which the 
reservoir conditions are determined from input data on the total mass flow and 
the stagnation enthalpy. The mean* stagnation enthalpy of the gas stream can be 
calculated from measurements of the electrical power input, tl\ermal losses to 
the facility cooling system, and mass flow. In this option, the reservoir 
temperature and pressure are calculated by a double iteration to match the input 
mass flow and enthalpy values, assuming equilibrium flow from the reservoir to 
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the throat. If the boundary layer is included, a correction is made for the ef- 
fect of the displacement thickness on the effective throat area. 


1.6 Boundary Layer 

The boundary layer in an arc heated wind tunnel can be either laminar or 
turbulent. According to a correlation of available boundary layer transition 
data, the layer is expected to remain laminar throughout the nozzle for most 
operating conditions of the existing NASA/ JSC arc heaters. NATA contains an ap- 
proximate laminar boundary layer calculation based upon an integral method de- 
vised by Cohen and Reshotko. The calculation includes the effects of the 
streamwise pressure gradient in the nozzle. The Cohen-Reshotko method utilizes 
analytical curvefits to relations among non-dimensional boundary layer parameters, 
based on similar boundary layer solutions. The curvefits employed in NATA in- 
clude the dependence of the parameters upon the Prandtl number, the viscosity- 
temperature index, and the hypersonic parameter. The boundary layer is assumed 
to start at a specified position upstream of the nozzle throat, and is computed 
step by step along with the inviscid flow solution. Beyond the throat, the in- 
viscid flow is coupled with the boundary layer through the effect of the displace- 
ment thickness on the effective area ratio. The coupled flow is stabilized with 
the aid of a computational artifice. In the case of flow in a rectangular chan- 
nel, two separate boundary layer solutions are computed, one for each pair of 
channel faces. Besides determining the displacement effects on the inviscid 
flow, the boundary layer solution yields predictions of the heat flux and shear 
stress on the nozzle or channel wall. 


1.7 Flow Solutions 

NATA provides frozen, equilibrium, and non-equilibrium flow solutions. 

These are the same flow options offered by the original CAL program (ref. 1) and 
the methods used in generating the solutions are largely identical with those 
used in the CAL program. In the frozen solution, the species mole fractions are 
held constant at their reservoir values. This type of solution approximates the 
actual non-equilibrium flow fairly well in cases with low reservoir pressure. 

In the equilibrium solution, the gas is assumed to be in a state of local thermo- 
chemical equilibrium at each point in the nozzle. The actual flow always departs 
radically from equilibrium in the supersonic region downstream of the throat. 
However, equilibrium flow is often a good approximation in the region upstream 
of the throat and in the throat region. 

The non-equilibrium solution is intended to model the actual flow as closely 
as possible within the basic approximation of quasi-one-dimensionality. The 
species concentrations are assumed to be governed by chemical rate equations. 

The basic method of solution is numerical integration of these rate equations to- 
gether with the differential equations formulating conservation of mass, momentum 
and energy. However, it is found that the system of difference equations is sta- 
ble only for extremely small step sizes when the flow is close to equilibrium. 
Because the flow starts from a state of equilibrium in the reservoir, the numeri- 
cal integration technique therefore cannot be used initially. Instead, the 
solution is started by treating the non-equilibrium flow as a perturbed equili- 
brium flow. The perturbation method is used until the departure from equilibrium 
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has become large enough to allow use of the numerical integration. Another ■ 
difficulty is that the sonic mass flux: for the non-equilibrium solution is not 
known until the solution has reached the sonic point. This problem is cir- 
cumvented by using an inverse method upstream of the throat and for a short dis- 
tance beyond the throat. 

When a two- temperature electronic non-equilibrium model is used, the electron 
temperature is one of the integration variables. In such cases, radiative energy 
losses and energy transfer between the electrons and the heavy particles in the 
gas are taken into account, and some of the electronic excited states may be 
treated as separate species governed by rate equations and not necessarily in 
equilibrium with the ground state. 


1.8 Model Conditions 

A major objective of NATA is to provide calculations of test conditions on 
models for comparison with experimental heat flux and pressure measurements. 

The code provides model calculations of two types: stagnation point conditions 

and conditions on blunt wedges. 


1.8.1 Stagnation Point 

The calculations of stagnation point model conditions begin with a normal 
shock solution. The gas immediately behind the shock can be assumed to be fro- 
zen (with mole fractions equal to those in the free stream ahead of the shock) 
or in chemical equilibrium. Either or both of these types of solution can be 
obtained in a given NATA run, regardless of whether the free stream solution is 
frozen, equilibrium, or non-equilibrium. The conditions at a point on the stag- 
nation streamline just outside the boundary layer on the model are then calcu- 
lated. The pressure at this point is the stagnation pressure, which can be 
compared with Pitot measurements. The stagnation point heat flux is then calcu- 
lated using a modification of the Fay-Riddell formula. Heat flux calculations 
are done for both hemispherical and flat-faced models, and for both an equilib- 
rium and a frozen boundary layer. Effects of surface catalytic efficiency can 
be included. The shock standoff distance is also calculated approximately. 

These calculations of stagnation conditions on models are done independently for 
the frozen shock and the equilibrium shock. 

Calculations of stagnation point conditions are normally done for axisym- 
metric models. However, the calculations can be done, instead, for two-dimen- 
sional models such as a cylinder with its axis normal to the direction of flow. 


1.8.2 Wedge Models 

The pressure and heat flux distributions on the flat surface of a blunt 
wedge are calculated using a modification of the results of the Cheng-Kemp theory. 
Effects of bluntness and the boundary layer displacement thickness are taken into 
account. For a given model location in a particular NATA solution, wedge condi- 
tions can be calculated for a series of input-specified leading edge radii and 
wedge angles of attack. 
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1.9 Supplementary Program Functions 


NATA contains several features designed to aid the user during the develop- 
ment of new gas models and under other abnormal circumstances. These features 
are not normally exercised during runs to produce flow solutions based on 
standard gas models. 


1.9.1 Equilibrium Gas Properties 

The code contains an option to compute only the reservoir conditions, i.e., 
the equilibrium thermodynamic properties and mole fractions. This features al- 
lows the user to determine the equilibrium equation of state of a gas mixture of 
specified composition, based on a particular set of assumed species properties. 
Normally, the transport properties in the reservoir are also computed. However, 
the transport property calculations can be suppressed, if desired, to avoid in- 
put of the required collision cross sections. 


1.9.2 Species Thermal Properties 

Another option produces tables of the free energy, enthalpy, specific heat, 
and entropy as functions of temperature for each of the species in a gas model. 
This feature provides, for example, a convenient means for testing proposed 
thermo fits for the species, and allows direct comparison of the properties, as 
computed by NATA, with other tabulations such as the JANAF tables. 


1.9.3 Transport Cross Section Edits 

A third option produces an edit of the steps in the transport cross section 
calculation for all of the species in the gas model. This edit is a useful aid 
to setting up or modifying cross section inputs for a gas model. In addition, 
NATA can be made to produce a deck of punched cards containing the averaged 
collision cross sections for all pairs of species in a gas model. In this form, 
the data can be read and used by other computer programs. 


1.9.4 Tape Output 

NATA contains provisions for writing selected results of the flow, boundary 
layer, and model calculations on a binary tape for subsequent processing by 
other programs. An auxiliary program (NATA/PLOT) is available for producing cer- 
tain types of plots of such data using SD-4060 or similar equipment. 


1.9.5 Error Processing 

When a large, multicase NATA job is run, one or more of the cases may fail 
because of input errors, inadequacy of standard values for the control parameters, 
or possible previously undetected coding errors. NATA contains numerous provi- 
sions to aid the user in identifying the cause of trouble in case of code failure, 
and to facilitate the running of large jobs in spite of the failure of some of 
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the cases. First, the program contains many validity checks designed to detect 
errors before they cause the program to execute an operation which is not al- 
lowed by the computer system, such as dividing by zero or calling a Fortran func- 
tion with an illegal argument. Such an operation would lead to immediate 
termination of the entire run. Second, when an error condition is detected, in 
most cases NATA prints a diagnostic message to identify its nature. Third, 
execution of the current case in the job is terminated. Fourth, a special sub- 
routine is called to print out most of the data stored in common, together with 
alphanumeric identifiers. The name of the subroutine in which the error was 
detected is also printed. Finally, the program proceeds to read the input data 
for the next case, and execution of the job continues. 


1.10 Limitations 

The calculations performed by NATA involve many approximations which limit 
the code's accuracy and applicability. The basic approximation in the inviscid 
flow solutions is that of quasi-one-dimensionality. Each flow variable is as- 
sumed to vary with the axial coordinate of the nozzle, but to be constant over 
any nozzle cross section. Actual flows in arc heated wind tunnels always show 
some radial non-uniformity. In some facilities, the non-uniformity is suffi- 
ciently moderate that the quasi-one-dimensional description can be considered to 
be a roughly valid and useful idealization. In others, the non-uniformity is 
so severe that the usefulness of an analysis based on radially uniform flow ap- 
pears doubtful. 

The boundary layer calculations involve many approximations related to the 
boundary layer starting condition, the method of solution, and the stabilization 
of the coupled boundary layer/inviscid flow problem beyond the throat. Also of 
course, the laminar boundary layer calculation is applicable only when the actual 
boundary layer has not undergone transition. 

It would be extremely difficult to carry through an error analysis for the 
code by evaluating and combining the effects of all of the approximations used. 
Such a calculation has not been attempted. Instead, results from the code have 
been compared with experimental data from arc heater facilities in which the 
radial non-uniformity is not very severe. It has been found that NATA predic- 
tions of static pressure and stagnation pressure generally agree with the 
measurements to within about 20 percent. Based on a limited number of compari- 
sons, it appears that NATA predictions of heat flux to a channel wall or a wedge 
model are too low by about 20 to 30 percent. NATA results for stagnation point 
heat transfer are roughly in agreement with experimental data for hemispherical 
models, when allowance is made for surface catalytic efficiency and for low- 
density effects which are not treated by the code. 


1.11 Organization of Report 

The remaining sections of this report document the theoretical relations 
and mathematical analyses upon which the NATA code is based. Section 2 treats 
the description of chemical species and chemical reactions. Section 3 deals 
with the calculation of gas transport properties. Section 4 explains how the 
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geometry of the nozzle or channel is described, and how the boundary layer dis- 
placement effect on the inviscid flow is formulated. Section 5 gives the 
analytical basis of the boundary layer calculation. Section 6 discusses the 
equilibrium and frozen inviscid flow calculations, and Section 7 the non-equilib- 
rium flow. Finally, Section 8 treats the calculations of conditions on models. 
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2. THERMODYNAMIC AND KINETIC MODELS 


The treatment of chemical equilibrium and reaction kinetics in NATA is still 
essentially the same as that in the Cornell Aeronautical Laboratory program 
(ref. 1). Section 2.1 explains how the code describes chemical species and 
formulates the conservation of chemical elements and electric charge. Section 2.2 
discusses the calculation of thermal properties of species and Section 2.3 the 
specification of chemical reaction rates. Simple models for vibrational non- 
equilibrium and for imperfections in the gas equation of state are described in 
Sections 2.4 and 2.5. 


2. 1 Elements and Chemical Species 

The gas stream whose flow is computed by the NATA code can consist of up to 
20 chemical species. The types of species- normally treated are atoms (including, 
in some instances, atoms in specific electronic excited states), atomic ions, 
diatomic molecules and molecular ions, linear triatomic molecules or ions (such 
as CO 2 ) , and the electron. The number of species included in a particular prob- 
lem will be denoted by n . The gas composition is expressed, in different parts 
of the calculation, in terms of either the mole fractions Xj or the molar concen- 
trations yj in units of moles of species j per gram of mixture. These quantities 
are related by 

Xj = W yj (!) 

in which W denotes the local mean molecular weight of the gas mixture in grams 
per mole. 

The number of chemical elements present in the gas is denoted by c. The 
chemical formulas of the species are represented by a matrix a, whose general 
element a- is the number of atoms of the j th element per molecule of the i th 
species. Thus, if the species is denoted by Mj and the element by E- , 
the chemical formula of the species may be represented by the equation 


M i ' X) ai > E i 

j = 1 


( 2 ) 


If ion species are present, the electron is included among the elements. In that 
case, positive ions are represented as compounds containing a negative number of 
electrons; for example, N 2 + is considered to be a compound N 2 e_i. With this con- 
vention, conservation of electric charge during reactions becomes a special case 
of conservation of the chemical elements. 

NATA calculations of thermochemical equilibrium* are carried out using a 
technique in which the concentrations of the "dependent species" are expressed 
in terms of those of the "independent species" (or "components"), where the 

*See Section 6. 
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number of independent species is equal to the number of elements. The advantage 
of this technique is that it reduces the number of equations, which must be solved 
simultaneously to determine the equilibrium species mole fractions, from n toe, 

For example, the high temperature model for air used in NATA contains 11 species 
but only 3 elements (n, O , and e ) , so that the number of equations is reduced 
from 11 to 3. 

The independent species or components must be chosen such that they are 
linearly independent combinations of the chemical elements, and must be placed 
at the beginning of the list of species for each gas model. To avoid computational 
problems, they should be species which are stable at low temperatures (i.e., at 
large expansion ratios downstream of the nozzle throat). However, if ion species 
are included in the model, the electron should be one of the components, and should 
be listed as species number 1. The ion species should all be placed at the end of 
the list of species. This arrangement allows the code to drop the charged parti- 
cles from the model for equilibrium flow when the equilibrium electron mole frac- 
tion becomes negligible (< 10 “ 20 ) . As an illustration of an arrangement satisfying 
these requirements, in the NATA models for air the components are e, N 2 » and 02 > 
in that order. 

Because the components head the list of species, the relation between them 
and the elements is specified by the square submatrix ajj of a ; - with i = 1 to c . 
Since the components are chosen to be linearly independent combinations of the 
elements, this submatrix is non-singular and has a unique inverse which will be 
denoted by • The inverse is defined by 

c 

y A ki «ij = s kj (3) 

i = 1 


where the Kronecker symbol <5^- is 1 for k = j and 0 otherwise. With the aid of 
the inverse matrix A^, the system of equations (2) for i = 1 to c can be solved 
to obtain the elements as linear combinations of the components: 


2 A ki M i 
i = 1 


yi a w E j - y] ^ ^ % 

i = 1 1=1 j = 1 


(4) 


Equation (4) can now be substituted into the remaining equations (2) for i = c + 1, 
c + 2 , . . . , n , to obtain expressions for the dependent species in terms of the 
components: 


M: 


y "i-c,k M k <* - c + l, . . . ,n) 
k = 1 


( 5 ) 
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where 


v i- c, k 


The matrix v gives the composition of the dependent species in terms of the com- 
ponents . 

The overall elemental composition of the gas mixture is invariant through 
all reactions, equilibrium or non-equilibrium. It is specified on input in terms 
of the mole fractions Xj c of the species making up the cold gas which is injected 
into the arc heater. These mole fractions and the chemical formulas for the cold 
species are used to calculate the number of gram-atoms of each element k per mole 
of the cold gas, denoted by Q k : 


Qk = X x j c “ik c (7) 

i = 1 


E 

j = 1 


a ij A jk 


(i = c + n) 


( 6 ) 


where aj k c is the number of atoms of the k tfl element in a molecule of the j th cold 
species, and n c is the number of cold species. 

For use in the thermochemical equilibrium calculations, the elemental com- 
position of the gas must be re-expressed in terms of the independent species. 
Because the elements can be expressed as linear combinations of the independent 
species by equation (4) , the elemental composition of the gas is equivalent to a 
composition in terms of the components. Let qj c denote the number of molecules 
of the j th component per molecule of the cold gas when the gas is considered to 
consist of the independent species. These quantities qj c are related to the Q k 
by ! 


Qk = X q > C a ' k (k = l ’ • • ’ > c) 

j = l 


( 8 ) 


Multiplication of this equation by A ki and summation over k gives 


q i ~ A ki Qk 1, • • • , c) 

k = 1 


( 9 ) 
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NATA computes the q? from equation (9) , and then normalizes them 


9f 


9i = 


j =1 


(i = 1, . . c) 


( 10 ) 


to obtain a set of composition coefficients whose sum is unity. If the gas 
actually consisted of the components and no other species, the q. would be the 
mole fractions of the components. 


The above definitions and relations can be used to derive a formula -for the 
mole fractions of the independent species in terms of the qj and the mole frac- 
tions of the (n-c) dependent species. This formula is used in the equilibrium 
calculations. Let Nj represent the number of molecules of species j in a gas 

n 

sample containing N = Y', Nj molecules altogether. Also, let Nj denote the 

number of molecules of the j th independent species that the sample would contain 
if all of the atoms making up molecules of the dependent species were rearranged 
to form molecules of the independent species. From the definition (5) of the 
*i-c,k matrix, 


N: 


N i 


N; V 


i-c,J 


<1=1 c) 


■ c + 1 


( 11 ) 


where the first term is the number of molecules of species j actually present and 
the sum represents the extra species j molecules that could be formed from the 
dependent species. From the definition of q. , 


9i 


N; 


■S N k * 


k = 1 

From (11) and (12), 


q i 2 [ Nk + 5^ N i *'i-c,k] = N j + 2 Ni 


k = 1 


i = c + 1 


i = c + 1 


( 12 ) 


(13) 


This equation may now be divided through by N to convert the subscripted N*s to 
mole fractions, X. It is noted that 
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£ 


x k = 


*-E 


x i 


(14) 


Also, 


the singly subscripted array 


is defined by 


* 



c 


c 



k = lr 


(15) 


Equation (13) then becomes 
n 

Xj = qj — X i f v i-c, j — ( v i- c — 1)-^ (16) 

i = c + 1 


which gives the mole fractions Xj of the independent species (j = 1 to c ) in terms 
of the <lj and the mole fractions 'X. ( i = c + 1 o') for the dependent species 


2.2 Thermodynamic Properties of Species 

The thermochemical equilibrium calculations for the reservoir and the nozzle 
flow calculations require data on the following properties for each chemical 
species (j ) in the gas: the molar enthalpy Hj , the chemical potential /x° and 

molar entropy S:° at a standard pressure p° , and the molar heat capacity c Pi • For 
a mixture of ideal gases, these thermal properties are functions only of temper- 
ature. 

The chemical potential may be defined (ref. 2, p. 283) as the partial deriva- 
tive of the free energy of a gas mixture with respect to the quantity of one of 
its constituents: 

/Xj = (dF/dNj) S; v>I j k (17) 

where the subscripts indicate that the entropy S , volume V, and quantities of the 
other constituents are held constant. The units of /xj depend upon those for F 
and Nj . The quantity measure Nj can be molecules, moles, grams, etc. On a molar 
basis, /xj can be calculated using the relation (ref. 2, p. 953) 

, / f i\ ° (18) 

Py = - R 0 T ^ + H j0 
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in which R 0 is the universal gas constant, fj is the molecular partition function 
for speciesij , defined as a sum over all energy states E] 



i 


and Hj*o is the enthalpy of formation of the species at standard conditions 
(usually zero temperature) . For an ideal gas 


(2 Trmj e 

n ; 


W <T> 


( 20 ) 


where f; nt is the partition function for the internal degrees of freedom. In the 
first factor on the right, which is the translational partition function, nj de- 
notes the number of particles of species j per unit volume and mj the particle 
mass; k is Boltzmann's constant, h Planck's constant, and e the base of natural 
logarithms. The partition function depends upon the pressure through nj . For 
convenience in calculations, equation (18) is rewritten in the form 


Mj = Mj° + R 0 T fn Pj (21) 

where pj denotes the partial pressure in atmospheres and pj° is the chemical po- 
tential at a partial pressure of 1 atmosphere (a function only of temperature) . 
From equations (18), (20), (21) and the ideal gas law 


Pj = n j k T 


it is possible to show that 
Mj j° = - R 0 T In 


'fj kT 


+ H; 


jO 


( 22 ) 


where p° = 1.01326 x 10^ dyne/cm^ (a pressure of 1 atmosphere expressed in 
absolute cgs units) and rj is a function only of the temperature: 


rj ^ njfj/e 


Because the molar internal energy of a species is given by (ref. 2, p. 335) 


'din f; 


Ej = R 0 T 


dr 


(23) 


(24) 


V, N; 
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one can show, by differentiating (22) , that the molar enthalpy H- = E- + R 0 T is 
given by 


H ) ~ H j0 _ _d_ 

R 0 T ~ T dT 


( o o 

fj - H j0 

R 0 T 


(25) 


° 

in which Hjq is the enthalpy of formation of the species. The molar enthalpy of 
the gas mixture is simply the sum of the species enthalpies weighed by the mole 
fractions : 



(26) 


The other required 
H: and [±- . The entropy 
287) ’ 


species thermodynamic properties can be calculated from 
Sj of the j th species is given by (ref. 2, p. 261 and 


e _ n □_ 

j T 

With (21) , this can be rewritten 

s i = s j - R ° k* Pj 

where pj is the species partial pressure in atmospheres, and 

, "i-4 

s j T 


(27) 


(28) 


(29) 


is the species entropy at the standard pressure (1 atm), a function only of 
temperature. The molar entropy of the gas mixture is given in terms of the 
species entropies by (ref. 2, p. 618) 


s . £ x jSj 

i - 1 


If p denotes the total pressure in atmospheres, (30a) can be rewritten, with the 
aid of the relation pj = Xjp between the partial pressures and mole fractions, in 
the form 
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n 

n 


7! x j s j° - R 0 P - % 

X) x i *" x i 

(30b) 

j = i 

j = 1 



in which the last term is the entropy of mixing. 

The molar heat capacity of the j th species is defined by 
d Hj 

Si - TF (31) 

The heat capacity of the gas mixture is then 
n 

c P - 2 si «» 

j = i 


Equations (25), (29) and (31) show that H- , Sj° and Cp- can all be calculated 
from the chemical potentials at the standard pressure. Wo methods are avail- 
able in the NATA code for calculating the chemical potentials themselves. Both 
are based fundamentally upon the statistical mechanical expression (18) for ftj in 
terms of the molecular partition function fj . One method, called the thermo-fit 
technique, relies upon accurate calculations of /i° by other computer programs. 

The data on as a function of temperature are curvefitted in the form 


o It o 

Fj HjO 1 pi 2 

- a j( l-^T) - bjT - - CjT 2 - - d jT 3 - 


R 0 T 

Then from (25) , 


— e.T 4 - k: 
4 J ’ 


Hj - Hj 0 
R 0 T 


+ bjT + cjT 2 + djT 3 + ejT" 


(33) 


(34) 


Thus, the parameters aj , bj , etc., are simply the coefficients of a fourth-degree 
polynomial curvefit to (Hj - H°q )/RqT. The thermo-fit method has the following 
advantages : 

1. The thermal property data generated by a specialized thermal-property 
computer program can be highly accurate, because the program can perform 
elaborate calculations, taking account of effects such as vibrational 
anharmonlcity and vibration-rotation coupling in molecules. 

2. The computer time required for evaluation of the curvefit formulas (33) 

and (34) is small. 
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3. The method is applicable to molecules which are too complicated to be 
treated by the other technique available in NATA, e.g., non-linear tri- 
atomic molecules and molecules containing four or more atoms. 

The curvefit coefficients aj , b- , etc., are determined by least-squares 
fitting of equation (34) to computed enthalpy values over a certain temperature 
range. Within the range of the data, the fit has good accuracy, but abo.e and 
below this range the accuracy can become very poor. If the range of the data is 
simply extended, the overall accuracy deteriorates because the polynomial form 
(34) is not sufficiently flexible. 

This limitation of the thermo-fit method can be circumvented by using the 
second NATA technique for calculating species thermal properties to obtain the 
properties at relatively low temperatures. The second technique, called the 
"physical model"*, is based upon approximate evaluation of the internal partition 
function fj nc in (20) for each species. For a monatomic species, the only in- 
ternal energy states are those of electronic excitation, and 


mt 



-E: /kT 

*i e 


(35) 


In this equation, Ej is the energy of the i th electronic state relative to the 
ground state (i = 1 ), and gj is the state degeneracy. In principle, the upper 
limit of summation L should be at the highest bound state below the effective 
ionization potential. In practice, results of good accuracy can be obtained by 
including only a few states of the lowest energy. High-lying states of. nearly 
equal energy can be lumped together by using an average Ej and summing the 
degeneracies. Array dimensions in NATA allow use of up to 10 states, including 
the ground state. 

For species containing more than one atom, the internal partition function 
sum includes states with various values of the vibrational-rotational energy. 

In the NATA physical model, the partition function for molecules is approximated 
as a product 

^int = ^rot ' ^vibr ‘ ^elec 

of rotational, vibrational, and electronic factors. This approximation neglects 
the variation of the molecular moment of inertia with the vibrational quantum 
number as well as the difference in vibrational frequencies between different 
electronic states of the molecule. As a further approximation, the vibrational 
energy states for each normal mode are approximated by those of a harmonic 
oscillator. These approximations to the partition function for molecules be- 
come inaccurate at high temperatures, but are expected to give good results at 
moderate temperatures where most of the molecules are in their electronic ground 
state and the degree of vibrational excitation is not too high. 


*lt is called the "harmonic-oscillator model" in reference 1 . 
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For linear molecules (including all diatomic species) , the rotational 
partition function, f rot , neglecting quantum effects as is appropriate for the 
temperatures of interest in NATA, is given by (ref. 2, p. 345 and 454) 

8rr 2 IkT kT 

W = T- - — ( 

a h x a h c Bq 

where I denotes the molecular moment of inertia about an axis, passing through 
the center of mass, perpendicular to the line joining the constituent atoms. 
The symmetry factor a is unity for asymmetric molecules and equal to 2 for 
symmetric molecules (with identical nuclei at opposite ends of the molecule) . 
The quantity Bq in the final expression of (37) is the rotational constant of 
the molecule in the ground vibrational state. For diatomic molecules, it is 
given in terns of tabulated spectroscopic constants by 


B 0 = B e - 



(38) 


where B e denotes the rotational constant for the equilibrium internuclear separa- 
tion and a e is the coefficient giving the dependence of the rotational constant 
upon the vibrational quantum number . The term a e is normally quite small and 

can be neglected without serious error within the context of the approximations 
being used. In terms of B 0 , the moment of inertia is given by 


I = h (39) 

8 v 2 c Bq 

Thus, (37) can be written 


f iot ^ 


(40a) 


where the characteristic rotational temperature d T is defined as 


hccr 

6 C = — j— B 0 = 1.43879 o Bq 


(40b) 


with Bq given in cm 


For non-linear triatomic molecules* (ref. 2, p. 454), 

1 /? 9 3/2 , 1/2 
ir 1 ' 2 (8 it 2 kT) (I A I B ^ 


‘rot 


a h3 


( 41 ) 


where I A , I 0 , and I c are the three principal moments of inertia. 


•The non-linear case is not treated in the present version of NATA. The formulas for this case are given with a view toward possible 
future development of the code. 
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For diatomic molecules, the vibrational partition function in the harmonic- 
oscillator approximation is (ref. 2, p. 345) 

«*,- (>- <- h -‘ /tT r4-' _9v/T r < 42 > 

where hvj is the energy difference between the ground state and the first vibra- 
tionally excited state, and 0 v the characteristic vibrational temperature 

0 V = ht/j/k (43) 


The vibrational term of a diatomic molecule is represented in the form (ref. 3, 
p. 92) 


G(v) = <u e 




y 


e 




(44) 


The vibrational constants <u e , ai e x e , and <u g y e have been determined from spectro- 
scopic studies for many diatomic molecules (ref. 3, p. 501-581). The vibrational 
excitation energy hvj is given by hv! = he [ G (1) — G(0) ] . Thus, from (43) and 

(44) 


e 


V 


he 

k 



- 2co e x e 




13 

4 


“eVe 


(45) 


It is clear from (44) that G (v) - G(v-l) varies with the vibrational quantum 
number v . For this reason, the present harmonic oscillator approximation (which 
neglects this variation) becomes inaccurate at high temperatures where several 
vibrational levels are populated to a significant degree. 

A triatomic molecule has three vibrational normal modes. The vibrational 
partition function for a non-linear triatomic species is 

<*■ n (>->T <“ 

k = 1 


In a linear triatomic molecule, there are only two rotational degrees of freedom 
and the "bending" vibrational mode is degenerate. In this case, the vibrational 
partition function is given by 

u,- n (4?: 

k = 1 


in which two of the 0^ , say $ v 2 = 0 V 3 , are equal to the vibrational temperature 

for the bending mode. 



The chemical potential at standard pressure can now be evaluated by sub- 
stituting the expressions for the internal partition function, based on the physi- 
cal model, into the general relations (20), (22) and (23). For monatomic species, 
fj nt is given by (35) , so that 


f'j’-HjO \ 3 /2nm k\ / k \ ' ? 5 

* u [—)* u br* 


(48) 


For molecular species, from (36),, 


Vj - Hj 0 
R 0 T 


2 TT m: k 


in 


+ in [ + — in T + in ( L t ) 


(49) 


+ in <f v ibr ) + In ( f elec^ j 

The chemical potential formulas for monatomic and diatomic species can be com- 
bined into a single equation involving the number nj of atoms in a molecule of 
the species : 


o o 1 

P) - H i0 _ _) 
R 0T = | 


where 


5 + 2(nj — 1) _ 0 . /t 

b; + in T - ( n: — 1 ) in ( 1 - e V J ) 

) 2 > 


in 


L i 


i = 1 


.. - E ii /kT 


*ij e 


(50) 


b; 


1 in ( 2rrm,k \ + tn(—\ - (n: - 1) in d- 


(51a) 


In these formulas, equations (35), (39) and (42) have been used to evaluate the 
electronic, rotational, and vibrational factors in the partition function. If 
numerical values are substituted for the physical constants, (51) becomes 


b ; = - 3-66505 + - in W: - (n, - 1) 6 ■ 

I 2 > ’ 1 


(51b) 


in which W- is the species molecular weight in grams per mole. 
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In the case of a linear triatomic species, from (49), (35), (39), and (47) 

4 


o o 
f*j ~ HjO 

R 0 T 


B; + 4- &» T — (1-e 
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'j ' 7 


i =1 


+ En 


8ii e 


-Ejj/kT 


i = 1 



(52) 


where 


i /2im:k \ / u \ 

B: = - l n [ — + £„ ( — „ e„ fl. 

’ 2 V h J / W " 


(53) 


The enthalpy can be calculated for monatomic and diatomic species by apply- 
ing equation (25) to (50) : 

H: - Hjfj 5 + 2(n:-l) 


~ H j Q 
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In the case of linear triatomic species, application of (25) to (52) gives 


Hj ~ H j0 7 

RqT ~ 2 



"vji 1 
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(55) 
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The molar heat capacity for monatomic and diatomic species is, from (31) 
and (54) , 


C P) 5 + 2 ( nj 1 ) 

= + (n* — 1) 

R 0 2 » 

S 1 s 3 - S 2 2 


m 2 ' 


«,/T 


«V T -.> 2 


where 
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i = 1 
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(56a) 


(56b) 


(56c) 
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For linear triatomic species, from (31) and (55), 
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(57) 


The thermal properties of a species can be specified using the thermo-fit, 
the physical model, or both techniques. In the case of a non-linear triatomic 
species or a species whose molecules contain four or more atoms, the programming 
of NATA does not provide a treatment based on the physical model. The properties 
of such species must be described using the thermo-fit at all temperatures. On 
the other hand, for monatomic species, the physical model is. accurate at all 
temperatures. Accordingly, only this model is used for the atoms and atomic ions 
whose properties are compiled into NATA. For most molecular species, both thermo- 
fit data and physical model properties are provided in the code. The physical 
model is used for temperatures up to a "switching" temperature T gw which is pre- 
set to 5000° K; and the thermo-fit is used for temperatures above T sw . To avoid 
disturbance of the flow solution by a small discontinuity in thermal properties 
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at the switching temperature, the transition between the two models is spread 
out over a temperature range 1000° K wide, between T sw - 500° K and T sw + 500° K. 
Within this range, the species thermal properties are calculated by mixing the 
results from the thermo-fit and the physical model; e.g., for the enthalpy, 

Hj = w TF • (Hj) Tp + w pM • (Hj) pw (58a) 

where 

(58b) 
(58c) 

and where the subscript TF denotes the thermo-fit and PM the physical model. 


TF 


V PM 


500 + (T-T sw ) 
1000 

500 - (T - T sw ) 
1000 


2.3 Reaction Rates 

A chemical reaction* involving some of the species present in the gas mix- 
ture can be denoted 


n 

X “*i M i 

k K 

t . 

x “i 

(59) 

j =1 

K ri 

i = i 



in which Mj represents a molecule of the j* species and the v- , v C. are 
stoichiometric coefficients. For the molecules which do not participate in the 
reaction as reactants, Vjj =0; for those which do not appear among the products, 
v\j - 0. The subscript i serves to distinguish the different reactions. The 
total number of reactions is denoted by r . 

The molar concentration of species j will be denoted by [Mj] . If (59) repre- 
sents an actual reaction mechanism, then the number of moles of forward reactions 
occurring per unit volume is given by 

n 

N fi = k fi J ~l 
k = 1 


[MiJ 


v ik 


(60a) 


‘See reference 4, Chapter XVII, for a background discussion of classical chemical kinetics. 
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in which the reaction rate constant kg is independent of the species concentra- 
tions but is, in general, a function of temperature. Correspondingly, for the 
reverse reaction 


n,i - n iMki 

k = i 


(60b) 


If (59) represents the net result of two or more steps in the actual reaction 
mechanism, then equations (60) are inapplicable. In such a case, if (60) were 
forced to fit experimental rate data, it would be found that the apparent rate 
constants k fi and k rf would vary with the species concentrations. 


From equations (59) and (60) , the net rate of change of the concentration 
[Mj ] of species j due to the i* reaction is 


d [Mj ] i 
dt 


= (vjj - vjj) kg [M k ] 


"Ik 


(61) 


k = 1 
n 


- (•'ij - v ij> k ri J"J [M k ] lk 


k = 1 


In equilibrium, the net rate of production is zero, so that, from (61), 


k ri = k fi J"! [M k ]( * 


(nk-nk) _ 

K: 


(62) 


k = 1 


where Kj denotes the equilibrium constant: 


Ki - 


(63) 


k = 1 
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The equilibrium constant for the reaction can be determined from, the condition 
(ref. 2, p. 953): 


n 



j = 1 



"ij ft 


(64) 


in which /ij denotes the chemical potential (21) for the j* species. Substitu- 
tion of (21) into (64) gives 



where 



( 66 ) 


For ideal gases, the partial pressures p- are related to the molar concentrations 
(Mj ] by 1 


Pj = [ Mj ] R 0 T 


(67) 


Hence, combination of (63) and (67) gives 


= 


(R 0 T) 


ft 


eX P 



where 


ft - ? ftj 


( 68 ) 


Because the chemical potentials at standard pressure (/tj° ) depend only on the 
temperature, the equilibrium constant K/ is a function only of temperature. 

The relation (62) between the forward and reverse rate constants for a re- 
action is termed "detailed balancing". Since k r j , kfj , and Kj depend only on 
the temperature, not the concentrations, this relation is valid even when the 
concentrations are out of equilibrium, as they normally are throughout a non- 
equilibrium flow solution. 
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Only one of the two rate constants for each reaction need be specified, as 
the other can be determined from the first using the detailed balancing relation 
(62) and the thermochemical expression (68) for the equilibrium constant. Con- 
ventionally, the forward rate is specified. In NATA, it is represented as a 
function of temperature by a curvefit of the form* 


m 


K fi = Aj 


10000° K 


-E ai /R 0 T 


(69) 


in which E a j is the activation energy for the reaction. Experimental data on 
rate constants are usually fitted using a form similar to this, but with T^i in 
place of (T/IOOOO)^. The specific form (69) has been adopted for use in NATA 
because, in reactions with large negative exponents qj, it permits fitting the 
experimental data with coefficients A] of smaller magnitude. When the conven- 
tional formula with T^i is used, some such reactions require Aj values larger 
than the limit (~1038) Q n floating point numbers in computers such as the UNIVAC 
1108. 


Many dissociation-recombination reactions involve a "third body", a particle 
which catalyzes the reaction by supplying part of the dissociation energy or 
carrying off part of the recombination energy. If the rate constants are the 
same for several reactions which differ only in the third-body species involved, 
the kinetic model can be simplified by combining all of these reactions into one. 
The third-body concentration in the combined reaction is equal to the sum of the 
concentrations of the actual third bodies. The coding of NATA allows this 
simplification. The third bodies are omitted from the v£j and matrices for 
the combined reaction but are listed separately. The rate-constant curvefit (69) 
for the actual three-body reactions is used. The effect of the combined concen- 
tration of all the third bodies for the reaction is taken into account by special 
coding. 


2. 4 Vibrational Non-equilibrium 

The calculations of species thermal properties performed by NATA (Section 2.2) 
normally assume that the vibrational degrees of freedom of molecules are in thermal 
equilibrium with the other degrees of freedom. However, the code also contains an 
option to calculate these properties on the assumption that the vibrational degrees 
of freedom are "frozen" at the reservoir temperature. This option, together with 
the normal, equilibrium property calculations, makes it possible to bracket the 
possible effects of vibrational non-equilibrium. When this option is used, the 
code calculates the species properties from the physical model at all temperatures. 
The thermo-fit cannot be used in this case because it is based on a curvefit to 
property calculations assuming complete equilibrium. 

The formulas for the species thermal properties in the case of frozen vibra- 
tion can be derived from the equations in Section 2.2 by assuming that the vibra- 
tional excitation temperature is constant and equal to the reservoir temperature, 
Tq. For diatomic molecules, the chemical potential (50) becomes 


‘Other forms are used for some of the reactions in the electronic non-equilibrium model for argon. 
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H-° - H-'q \ 5 + 2(n- - 1) " _q ./Xn\ 
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( 70 ) 


The enthalpy (54) becomes 
Hj — Hjq° 5 +2(nj — 1) 
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(71) 


Because the vibrational contribution to the enthalpy Hj is constant, the heat 
capacity C p j (31) does not contain a vibrational term. The species entropy at 
standard pressure is given by 


S:° 5 + 2(n: — 1) 

— — = b: + (1+M 
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For a linear triatomic species, the thermal properties in the case of frozen 
vibration are given by 
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The specific heat again lacks a vibrational contribution. 


2.5 Gas Imperfections 

The CAL program upon which NATA is based (ref. 1) contained an option to in- 
clude the effects of gas imperfections due to the finite volume of the molecules. 
Such effects are negligible under the conditions to which NATA is normally applied. 
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However, to provide for possible, presently unforeseen, future applications of 
NATA to nozzle flows with extremely high reservoir pressures*, this option has 
been retained. 

The CAL treatment of gas imperfections is based upon the van der Waals 
equation of state (ref. 2, p. 581) 


R 0T _ 

V-b 0 V? 


(76) 


in which V represents the molar volume 


V = 


w 

p 


(77) 


The quantity bo is the molar volume from which the centers of the molecules are 
excluded because of the finite size of the molecules. If the molecules are 
assumed to be spherical, then 


b 0 



(78) 


where Nq is Avogadro’s number and & denotes the molecular diameter. The second 
term in (76) represents the effects of attraction between the molecules. In NATA, 
the coefficient a in this term is assumed to be zero, because the attractive 
forces have negligible effects in the high temperature flows to which NATA is 
applied. The remaining term is written 

pR 0 T 

p = (79a) 

W 


where the effective density ~ is defined as 

P “ ~ (79b) 

^ b 0P 

w° 

The molecular weight W° in (79b) is written with a superscript o to signify that 
the coefficient b 0 /w° is read as an input, and is assumed to be constant through- 
out the solution. The molecular weight W in (79a), however, is the local value, 
which varies in equilibrium and non-equilibrium solutions. 

The van der Waals equation (76) is inaccurate at high densities where V be- 
comes comparable with bo . Thus, equations (79) provide a valid approximation to 
the effect of finite molecular volume only so long as the term bop/W° is small 
compared with unity, say, less than about 0.1. 


•In air, these effects are of the order of 3 percent when the pressure in atmospheres is about 0.1 times the reservoir temperature in 
degrees K. 
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Because molecular attractive forces are neglected in (79) , there is no 
correction to the ideal-gas internal energy. However, the enthalpy is given by 


H 


p 

E + — = E + 
P 


P_ 

P 


R 0 T 

w 


(80) 


In terms of the enthalpy H ;cleal for an ideal gas. 



(81) 


In NATA, the computation of reservoir conditions is carried out for speci- 
fied temperature and pressure. The composition is assumed to be the same as for 
an ideal gas mixture, but the enthalpy and density are calculated from (81) and 
(79). These relations are also used throughout the frozen and equilibrium flow 
solutions, and in the initial portion of the non-equilibrium solution which is 
calculated by perturbing the equilibrium solution. The correction for gas im- 
perfection is not used during the non-equilibrium integration because, in 
practice, the flow remains near equilibrium (and thus is calculated by the 
perturbation method) until the density has dropped to values at which the correc- 
tion is negligible. 


Section 5.3 of reference 1 can be consulted for further discussion of this 
option. For air, this reference recommends a value of a = 2.6 A = 2.6 x 10~® cm 
for use in equation (78). 
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3. TRANSPORT PROPERTIES 


NATA requires transport property values for use in the laminar boundary 
layer calculations and in calculations of heat transfer to models. Some trans- 
port properties and related quantities are also printed out. All of the trans- 
port coefficients required in NATA are computed internally by a subroutine 
package which was adapted from a separate, previously developed transport proper- 
ty program. The basic data required by the calculations are compiled into the 
code for the standard gas species. Such data can also be set in the input. 

To facilitate the use of the code in computing inviscid flow solutions for 
gas models involving non-standard species whose transport cross sections may be 
unknown, all transport property calculations can be suppressed by input of a 
single control variable (N0TRAN) . This procedure also suppresses the boundary 
layer calculation and calculations of stagnation point heat transfer and of 
conditions on wedge models. 

According to kinetic theory (ref. 5) the transport coefficients in a mix- 
ture of gases depend upon the scattering cross sections Oij (defined below) for 
collisions between pairs of molecules. The Chapman-Enskog theory provides ex- 
plicit, but extremely complicated, formulas for the viscosity, thermal conduc- 
tivity, binary diffusion coefficients, and other transport properties in terms 
of the collision cross sections, the particle masses, the species mole fractions, 
and the temperature. Thus, in principle, the problem of calculating the trans- 
port coefficients for a given mixture consists of two parts: first, the deter- 

mination of the collision cross sections Oij for all possible pairs of species 
(i, j ); and second, the evaluation of the Chapman-Enskog formulas. 

The collision cross sections can be determined or estimated in many dif- 
ferent ways. They can be measured directly in molecular-beam scattering experi- 
ments or can be determined indirectly from the analysis of data on various gas 
properties. For simple systems they can also be obtained from ab initio or 
semi-empirical quantum mechanical calculations. 

For collisions between heavy particles, it is generally advantageous to 
express the collision cross sections liij in terms of the interaction potential 
4 > ij between the particles. The interaction potential can then often be approxi- 
mated by a simple empirical form containing one or more adjustable parameters 
which may be chosen to fit the experimental data. For example, one type of 
model frequently used in this way is the Lennard-Jones (6-12) potential 

♦- “[(t ) 12 - (tT) 

Here, r denotes the center-to-center separation of the two molecules, e and a 
are the adjustable parameters, and 4> is the energy of interaction. The param- 
eters can be evaluated by fitting theoretical predictions based on the models to 
experimental data on some material property, for example, the viscosity, thermal 
conductivity, binary diffusion coefficient, or second virial coefficient. 
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Once the interaction potential <f>\] between a pair of atoms or molecules has 
been determined, the averaged collision cross sections Qij^ ’ required in the 
transport calculations can be computed by classical mechanics. In general this 
step requires a lengthy numerical computation; however, this computation has al- 
ready been carried out for a number of commonly used forms for the interaction 
potential, and tabulated values of the averaged cross sections as functions of 
temperature are available in the literature. 

In the case of a pair of atoms with unpaired valence electrons, quantum 
mechanical calculations show that there can be several different interaction po- 
tentials which occur with different probabilities. In this case the total cross 
section may be obtained simply by calculating the cross sections independently 
for each potential and then averaging the resulting cross sections with the ap- 
propriate probabilities to obtain the total cross section. 

For collisions involving electrons, classical mechanics is not applicable 
and quantum calculations are required to relate the collision cross sections flij 
to the interaction potentials. In this case it is often simpler to analyze the 
data directly in terms of the collision cross sections, without going through 
the intermediate step of determining the interaction potentials . 

After the cross sections s ^ have been determined for all pairs of 

species in the gas, the transport properties may be computed by evaluating 
several large determinants. This step is straightforward, but can require a 
considerable amount of computer time in the case of a gas mixture containing a 
large number of species. 


In NATA, the compiled-in transport property data consist of parameter 
values and tables for direct calculation of the averaged cross sections 
The interaction potentials themselves are not used. Thus, the step of calculating 
the cross sections from the interaction potentials is bypassed. In addition, the 
amount of computation required to evaluate the mixture transport properties is 
greatly reduced by using approximate formulas, developed by Yos (ref. 7), in 
place of the full formulas of the first Chapman-Enskog approximation. 


3.1 Basic Equations 

The properties generated by the transport 'property calculations in NATA are 
the mixture viscosity p , the electrical conductivity a , the frozen Prandtl 
number 

Nftf = (83) 

and the atom-molecule Lewis number 

_ P C P f ° am (84) 

Le K f 

where Cpf is the frozen specific heat of the mixture at constant pressure, Kf is 
the frozen thermal conductivity, p is the mixture density, and D am is the binary 
atom-molecule diffusion coefficient for the mixture. The calculations of these 
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properties are based on the mixture composition as given by the equilibrium or 
non-equilibrium flow solution, or by the equilibrium normal shock solution, 
as appropriate. 


The viscosity [l and the translational component K tr of the thermal conductiv- 
ity are computed from an approximation to the first Chapman-Enskog formulas 
(ref. 5), developed by Yos (ref. 7): 


Yj XiAAjW + a (ct) ) 
i= 1 

: (85) 

n 

1 _ 5(a) Y, Xj/(A>) + i (a) ) 
i= 1 


Here a may represent either the viscosity jl or the translational thermal con- 
ductivity K tr ; n is the number of species in the gas; and X; is the mole fraction 
of the i* species. Also, 
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For the viscosity, the quantities and A|j® appearing in equations (86) 

and (87) are defined by 
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( 88 ) 


while for the translational thermal conductivity they are defined by 
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In these equations, Wj is the molecular weight of the i th species, Nq = 6.0225 x 
1023 molecules/mole and k = 1.3805 x 10 - 16 ergs/°K are the usual unit conversion 
factors. 


A-.d) = - /_ L_L_ 
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(90) 


(91) 


and the Ojj® s ) are averaged collision cross sections for the collisions be- 
tween species i and j , which are supplied as input to_the transport property 
calculations. Several different definitions of the s ) symbols have 

appeared in the literature.* The one used here is 


i] 
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where 
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1(1 +0/(1 + 2 £-(- 1) £ )] j (i - cos £ y) 4 tt (7jj sin X ^ X 


Here a jj = ctjj > s) is the differential scattering cross section for the pair 
i - j , y is the scattering angle in the center of mass system, g is the relative 
velocity of the colliding particles, and y = [m;mj/ 2 (mi + mj)kTr' g is the 
reduced velocity. In NATA, the cross sections are calculated as func- 

tions of temperature and gas composition for each pair of species in the mixture 
from basic cross section data, which are either in tabular form or are given as 
simple analytical functions of temperature and composition. 

*The quantity designated here as is called rr <j 2 • s ^* in ref. 5 and n SI jd®. s ) in ref. 6. 
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The frozen thermal conductivity Kf appearing in (83) is calculated from the 
modified Eucken approximation (ref. 8): 

K f = K tr + K int < 94 > 

where K u is the translational component of the thermal conductivity given by 
equations (85) to (91) , 



is the component of the thermal conductivity resulting from the internal ex- 
citation energy of the molecules, and c pi is the specific heat at constant pres- 
sure for the i th species. The frozen specific heat c pf for the mixture is then 
given by 


V = 


1 
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(96) 


where 



(97) 


is the average molecular weight, and XjW-/W is the mass fraction of species i 
in the mixture. 


For the electrical conductivity, the NATA code uses an expansion of the 
first Chapman-Enskog approximation (ref. 5) to lowest order in the ratio (m e /m; ) 
between the electron and atom masses, 
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where e is the electronic charge 
that the term X e Ae e ^ is to be 
Chapman-Enskog approximation. 


and the 
omitted 


(98) 


prime on the summation sign indicates 
from the sum. The complete first 


D; 


P N 0 A ij 


( 1 ) 


(99) 


is used for the atom-molecule diffusion coefficient required in the Lewis number 
calculations. 
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3.2 Accuracy of the Calculations 


The accuracy of equations (85) to (91) , which are used for the calculation 
of the viscosity and translational thermal conductivity in NATA, has been dis- 
cussed in reference 7, and it is shown that these equations generally give re- 
sults within a percent or less of the first Chapman-Enskog approximation, with 
considerable savings in computation time. The accuracy of the first Chapman- 
Enskog approximation itself depends primarily on the temperature dependence of 
the collision cross sections s ) for the various species in the gas. The 

approximation is exact for cross sections which vary inversely as \ff (i.e., 
constant collision frequency, or force proportional to the inverse fifth power 
of the internuclear separation) , and becomes progressively poorer as the 
temperature dependence of the cross sections deviates further from this law. 

For neutral-neutral and neutral-ion collisions, where the cross section depen- 
dence is close to 1 VT", the first Chapman-Enskog approximation generally gives 
transport properties within a few percent of the exact solution (ref. 5). Be- 
cause the collision cross sections are known only within about 10 to 

20 percent at best, for the species considered in the NATA code, the additional 
error introduced into the calculated transport properties by the use of the first 
Chapman-Enskog approximation for neutral-neutral and neutral-ion collisions is 
negligible. 

The effects of Coulomb collisions, for which the cross sections are propor- 
tional to l/T 2 , are not given very accurately by the first Chapman-Enskog 
approximation; for example, in the limiting case of a fully ionized gas, the 
transport properties calculated from this approximation differ by about a factor 
of two from those obtained from an exact solution of the Boltzmann equation 
(ref. 9). Coulomb collisions are treated in NATA using effective cross sections 
which are chosen to make the transport properties calculated from equations (85) 
to (99) agree with the exact solutions for a fully ionized gas at low pressures 
(ref. 9). The transport properties calculated from equations (85) to (99) for 
partially ionized gases using these effective Coulomb cross sections are then 
found to agree with accurate calculations for partially ionized gases (refs. 10 
to 12) within about 10 to 20 percent at all degrees of ionization, which is again 
within the accuracy of the available cross section data. 

The effective cross sections used for Coulomb collisions in the NATA code 
are given by the following equations: 
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(100a) 


n l2 (1 ’ = 2 ' 71 Q c 

n 22 ( 1 ’ 1) = 12.7 Q c 
Bjj = 1.5625 


0 12 (2,2) = 5.44 Q c 
0 2 2 (2,2) = 21.76 Q c 


where the subscripts e , 1 
and doubly charged Ions, 


Q c = 



l n (f A) 


and 2 represent respectively electrons and singly 


(100b) 


e is the electronic charge, 
3 (kT) 3/2 

= T e 3 (7rne )l/2 


(100c) 


is the Debye shielding parameter introduced by Spitzer (ref. 9), n e is the num- 
ber of free electrons per unit volume, and f is a correction factor to take 
account of the reduced effectiveness of the shielding at high electron densities. 
The correct form for the factor f is not known at the present time and a number 
of different approximations have been suggested for it (refs. 13 to 16). For 
electron number densities below about lO^- 7 to 10l8/cm 3 f is and all approxi- 
mations for the transport properties (refs. 13 to 16) agree within about 10 to 
20 percent; at higher electron densities, however, substantial differences 
between the different approximations begin to appear, and the correct treatment 
of the problem has not yet been established. The NATA code uses the approxima- 
tion suggested by Yos (ref. 7): 


f 



2 

64 n e 1/2 

n 3 

9 kT e 


(lOOd) 


However, as noted above, the errors in this approximation may be large at the 
higher electron densities. 

Errors in the first Chapman-Enskog approximation may also be large for 
electron-neutral collisions, ranging up to about a factor of three in the case 
of electron-argon atom collisions (ref . 10) , because of the strong temperature 
dependence of the cross sections which results from the Ramsauer effect in argon. 
For this case also, the NATA code makes use of effective collision cross sec- 
tions in equations (83) to (99) to calculate the transport properties, with the 
cross sections being chosen to match the experimental electrical conductivity 
data as well as possible,* For those cases in which comparisons have been made, 
this approach yields transport properties which agree with the available experi- 
mental data within experimental error, which is generally of the order of 20 to 
30 percent; however, larger errors than this are of course possible in other 
cases. 


•Note that electrical conductivity is the only transport property which is significantly dependent on the electron-neutral cross sections. 
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The accuracy of the modified Eucken approximation used in the thermal con- 
ductivity calculations for the NATA code has been examined in detail by Mason 
and Monchick (ref. 17). Errors in this approximation may approach 10 to 20 per- 
cent in cases where there is a rapid exchange of energy between excitation and 
translational states through inelastic collisions, as is normally the case for 
rotational excitation in low temperature gases; however, the errors become 
smaller when the exchange is less rapid, and are negligible in cases where 20 or 
more collisions are required for the exchange of energy between internal excita- 
tion and translation. Because the latter situation is normally the case for high 
temperature gases, errors resulting from the use of the modified Eucken approxi- 
mation should generally be a few percent or less in the NATA code, though errors 
of as much as 10 to 20 percent may arise from this source in exceptional cases, 
where the inelastic collision cross sections for some of the important internal 
states of a species are large. For atomic gases, of course, internal excitation 
does not ordinarily contribute very significantly to the thermal conductivity 
(ref. 18) , so that in such cases the modified Eucken approximation is necessarily 
valid, regardless of the behavior of the inelastic cross sections. 

In summary, the accuracy of the transport property calculations in the NATA 
code appears to be governed in almost all cases by the accuracy of the data on 
atomic and molecular collision cross sections used in the calculations. For the 
cross section data compiled into the code, this accuracy should generally be in 
the range from about 20 to 40 percent, except in the case of carbon containing 
mixtures, where only rough cross section estimates are available for some species 
and errors may consequently range up to as much as a factor of two or more on 
occasion. 

The discussion of the accuracy of the transport property calculations in 
NATA has not included any errors in the calculated transport property values 
which may arise from errors in the gas composition data used as input to these 
calculations. When such errors are important, they will of course decrease the 
absolute accuracy of the calculated transport properties and may need to be con- 
sidered in evaluating the overall accuracy of the code results. 


3. 3 Method of Calculation 

The formulas in Section 3.1 give the desired transport properties in terms 
of the temperature T , the species mole fractions Xj , molecular weights Wj , speci- 
fic heats cpj , and the average collision cross sections Oij^’ (T) . All of 

these data, except the cross sections, are provided by the gas dynamic calcula- 
tions in NATA. 

The cross sections for all of the pairs of species are computed in a series 
of steps. First, the cross sections for all pairs are set to zero. Then, in 
each step, the values of ^ , iK 2, 2 ^ , andB*!^ 1, ^ are computed by a particu- 

lar method (or option) with a particular set of parameter values ,_ and these 
values are added to the corresponding cross sections , Ojj( 2 > 2 ) , and 

Bjfj 1) for each pair of species (i, j ) to which the step is applicable. 

The information concerning the applicability of steps to species pairs is stored 
in index arrays. If only one step of the cross section calculation is applicable 
to a particular species pair, then the cross sections for that pair are the 
values computed during that step. If more than one step is applicable to the 
pair, the cross sections for' the pair are built up by adding contributions from 
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the different steps. This procedure provides considerable flexibility in the 
representation of the cross sections. If the cross sections are poorly known 
for several minor pairs of species , but are considered likely to be roughly the 
same for all pairs, then the cross sections for all of these pairs can be set 
in a single step. 

The code contains default provisions for estimating some cross sections if 
they are not specified explicitly in the compiled-in data or the input. If none 
of the specified steps is applicable to a particular pair, and if both of the 
species are ions, then the effective Coulomb cross sections (100) are used. If 
one species is neutral and the other ionized, the formula 


0..«, s) = a (£, s) T — 0.4 


( 101 ) 


is the default option. The constants A®> are compiled into the code. If 
both species are neutral and are unlike (not the same species) , the cross sec- 
tions are estimated using the simple mixing rule 



s) 





( 102 ) 


However, if cross section data are not specified for like-like collisions of a 
neutral species, the code does not attempt to provide estimates of the cross 
sections, but returns an error message and terminates the case. 

The compiled-in data specify steps for calculating the cross sections for 
the like-like interactions of all of the standard species, and for those unlike 
interactions for which experimental or theoretical cross sections are available 
in the literature. None of the standard gas models requires all of these steps. 
For example, the air models do not need the steps specified for calculating the 
cross sections in helium. At the beginning of each NATA case, the code edits 
the steps of the cross section calculation to eliminate those which are not 
needed in the current gas model, and to insert required but unspecified steps 
in accordance with the previously described default options. A printed summary 
of the steps in the unedited and edited cross section calculations can be ob- 
tained by setting a control input. 


3.4 Cross Section Models 

NATA contains 12 methods or options for calculating 0(1. 1) t n^> ^ , and 
B in a particular step. Each of these options is briefly described below. 
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3.4.1 Coulomb Cross Sections 


In this option, the cross sections are 
Q (1 > = 0.8 C x Q c 

n (2 > 2) = C 2 Q c < 103 > 

B* = 1.5625 

where Q and C 2 are specified constants, and Q c is given by equations (100b) 
and ( 100 c) . 


3.4.2 Exponential Potential 

Here, the cross sections are obtained from Monchick's (ref. 19) collision 
integrals for the exponential potential, 4>=Ae~ T ^P , which are compiled into the 
code. The input parameters are A/k and p . 


3.4.3 Charge Exchange Cross Section 


In this option, ^ and B* are calculated for a resonant charge ex- 

change cross section of the form 


( 2 , 2 ) 


Q ex = (A - B log 10 v) 2 

where v is the relative velocity. 0 

input parameters are A , B , and the atomic weight. The ^ and B*j 

puted can either replace the value computed in earlier steps of the cross sec- 
tion calculation or be added to them. 


(104) 

is not calculated in this case. The 

com- 


3.4.4 Tabulated Cross Section 

Here the cross sections are given in tabular form as functions of tempera- 
ture. 


3.4.5 Power Law Potential 


This option calculates collision cross sections for an inverse power law 
potential, <f> = Ar _J ? , based on the analysis of Kihara, Taylor, and Hirschf elder 
(ref. 20). The required inputs are 7 and the quantities 


Oi = ttT ^3 — ^ r/ 2/3 A* 1 ) (rj) 


2/i) 


Aj = 0- - -0 A< 2 > ( t ) ) / A^> ( 7 ,) (105) 
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where (jy) and A ®(j/) are tabulated functions which are given for both 

attractive and repulsive potentials in reference 20. 


3.4.6 Lennard-Jones (6-12) Potential 

This option calculates the cross sections for the Lennard-Jones (6-12) po- 
tential, equation (82). The input parameters- are f/k and a . The tabulated 
collision integrals are compiled into the code. 


3.4.7 Transferred Cross Sections 


This option allows cross sections calculated for one pair of species to be 
used also for other pairs, possibly with a constant multiplying factor. The 
cross sections are calculated from the formulas 


Q* 1 . « = Cj O-jd. « 

a (2 ’ 2) = Cj c 2 2) (io6) 


where the C k are constant factors and the subscript ij indicates the pair for 
which the cross sections were previously calculated. 


3.4.8 Empirical Mixing Rule 

Here the cross sections for a pair of unlike species (i, j ) are calculated 
from the empirical mixing rule (102) . These calculated cross sections are then 
added to the previous values for the pair. 


3.4.9 Fairing Option 


This option modifies the previously calculated cross section values for a 
species pair according to the formula 


5 «- ■> 

new 


-(£, s) 

f(T)0 o id 


(107) 


where f(T) is a linear fairing factor given by 


f(T) = max 


0, min I 1, 


T-T n 


t 2 -t 0 


(108) 


This option allows the use of different forms for the cross section in different 
parts of the temperature range, with a smooth transition between them. 
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3.4.10 Generalized Mixing Rule 


This option is a generalization of the empirical mixing rule (102) in which 
the cross sections are calculated from the formula 






(109) 


where (i, j ) and (m, n) are any specified species pairs. 


3.4.11 Same-pair Transfer 

This option calculates one of the averaged collision cross sections 
for a pair of species from the previously calculated value of a different 
for the pair. In terms of the notation 


5*® = o^ 1 . » 

«ii (2) - % (2 ' 2) 


o..(3) = B *. a-* 1 * » 

1J 1 ) 


( 110 ) 


the new cross section is calculated from the formula 

%j (m) = c %j (n) (111) 

where m and n are two specified integers in the range from 1 to 3, and C is a 
constant. 


3.4.12 Multiplication by a Constant 


This option multiplies previously calculated values of the collision cross 
sections for a pair of species by a constant factor, according to the formulas 


Si,' 1 ' - c, »] oU 

3 - C, c 2 [Sj/ 2 ' 21 ] old (112) 


®ij ^3 ^®ij^old 


This is the same as the option described in Section 3.4.7, except that here the 
cross sections for a species pair are obtained from previously calculated values 
for the same pair instead of those for a different pair. 
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4. NOZZLE AND CHANNEL GEOMETRIES 


4.1 Uses of Geometric Information in NATA 

The NATA code uses information about nozzle and channel geometries for the 
following three distinct purposes. 

4.1.1 Area Ratio Calculations 

In the equilibrium and frozen flow solutions, the gas temperature is taken 
as the independent variable. The code decrements the temperature repeatedly, 
starting at the reservoir value, and at each point calculates all the other flow 
variables, including the effective area ratio: 

A e = A e /A e (113) 


Here, Ag is the effective cross sectional area of the inviscid flow, while A' 
is the cross sectional area of the flow at the sonic point. The nozzle geometry 
is then used to determine the position in the nozzle, as specified by the axial 
coordinate x , corresponding to this effective area ratio. If the boundary layer 
on the nozzle wall is neglected, this solution for x is straightforward. If 
the boundary layer is included, the solution must be based on a relation between 
the effective area ratio A e , the boundary layer displacement thickness 8 , and 
the geometric area ratio. 

A g = A g (x)/A ^0 (114) 

In (114) , is the geometric cross sectional area at position x , and A^q is the 
cross sectional area at the throat. 


The non-equilibrium flow solution is started just downstream of the reser- 
voir as a perturbed equilibrium solution. The gas temperature is again the in- 
dependent variable, and the geometric description of the nozzle is used to de- 
termine the axial coordinate corresponding to each computed flow condition. 

When the perturbations exceed a certain limit, the solution is continued by in- 
tegration of the non-equilibrium rate equations, using x as the independent 
variable. In the region of the non-equilibrium integration, the geometric 
description is used to calculate the geometric area ratio for given values of x . 
If the boundary layer is included in the solution, it is also necessary to 
calculate the effective area ratio from the geometric area ratio and the dis- 
placement thickness. 


4.1.2 Boundary Layer Calculations 

The overall geometry of the nozzle or channel wall determines the conver- 
gence or divergence of streamlines in the boundary layer, and thus affects the 
rate of buildup of boundary layer thickness. For example, if the streamlines 
diverge, the gas flowing in the boundary layer expands laterally, and as a re- 
sult the thickness increases more slowly than it would if the streamlines were 
parallel. 
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4.1.3 Specification of Model Points and Channel Points 

NATA can provide calculations of stagnation-point conditions on axisymmetrie 
or two-dimensional models, and of the conditions along the surface of a blunted 
wedge model. One way of specifying the locations in the flow at which the lead- 
ing edges of such models will be placed is to input a list of test section 
diameters. Similarly, points at which flow calculations are desired in a rec- 
tangular channel can be specified by inputting a .list of values of the channel 
width. In each such case, the nozzle or channel geometry must be used to deter- 
mine the values of x corresponding to the input diameters or widths. 


4.2 B asic Geometric Options 

NATA has been programmed to provide flow calculations for three types of 
nozzle geometry: 

(1) axisymmetrie nozzle 

(2) two-dimensional nozzle 

(3) rectangular channel. 

In all three cases, the nozzle shape is specified by means of curvefits to noz- 
zle wall profiles. Each such profile may be described by a function y (x) , where 
y is the perpendicular distance from the nozzle axis to the wall, and x is an 
axial coordinate, zero at the geometric throat and increasing in the downstream 
direction. 

In the case of an axisymmetrie nozzle, the geometric cross sectional area 
of the nozzle at a station x is given by 

Ag = n-[y(x)] 2 (115) 

Thus, the geometric area ratio is 


ryWr 

L yo J 


(116) 


where y 0 = y(0). For a two-dimensional nozzle, the cross-sectional area per unit 
length in the z -direction is 

Ag = 2 y(x) (117) 


and the geometric area ratio is 


( 118 ) 


Description of a rectangular channel requires 
The cross sectional area at axial position x is 


two profiles, y( x ) and z(x) 
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= 4 y(x) z (x) 


(119) 


and thus the geometric area ratio is 


y(x) 2 (x) 
?0 2 0 


( 120 ) 


It is assumed in equation (120), and in the code, that both the profiles for a 
channel have their. minima at the same axial location, x = 0 . 


4.3 Profile Description 

Each profile used in NATA is represented by an analytical curvefit contain- 
ing up to 12 sections. The profile must be continuous and must have a continu- 
ous first derivative. At least two sections must lie upstream of the throat, 
and at least two must lie downstream. The throat must be a section boundary. A 
NATA case involves either one or two profiles, depending upon whether the flow is 
in a nozzle or a channel. 


Each section in a profile fit can have one of three available forms: 

(1) Straight Line 
y (x) = Pj + P 2 x 

(2) Circular Arc Convex Downward 

y(x) = P X - /p 3 2 - (x-P 2 ) 2 

(3) Circular Arc Concave Downward 

y(x) = P : + 7 p 3 2 - (x-P 2 ) 2 

In the second and third forms, P3 is the radius of the circular arc andP 2 , Pi 
are the x and y coordinates, respectively, of the circle center. 

The parameter values Pj , P2 , P3 for each section must be chosen so that the 
entire profile fit represents the given nozzle or channel profile with adequate 
accuracy, subject to the following conditions: 

y(0) = R 0 


( 121 ) 


( 122 ) 


(123) 


dy/ dx 
dy/ dx 
dy/ dx 


continuous everywhere 


< 0 

for 

x < 0 

> 0 

for 

x > 0 
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where Ro is the throat radius in centimeters. The last three conditions imply 
that dy/ dx = 0 at x = 0 . A separate computer program called NOZFIT has been de- 
veloped for computing the parameters of NATA nozzle profile curvefits from the 
kinds of data available on engineering drawings, i.e. , lengths, angles, and 
radii of curvature. 

An additional constraint on the profile fits has been established by experi- 
ence: There must not be any region of non-zero length in which the geometric 

area ratio is constant, or nearly constant, and approximately equal to 1. If 
such a region is present in the nozzle fit, the non-equilibrium solution usually 
becomes unstable in that region. Thus, if the actual nozzle being represented 
has a finite region of constant area (i.e., a straight tube) at the throat, it is 
necessary to modify the curvefit slightly in order to obtain reliable operation 
of the non-equilibrium solution. In practice, representing such a region as a 
converging sectioh with a convergence angle of 3 degrees has given satisfactory 
results in most cases. 


4.4 Relations Between the Geometric and Effective Area Ratios 


The laminar boundary layer on the nozzle wall displaces the streamlines in 
the inviscid flow, thus changing the effective nozzle geometry. The amount by 
which the effective boundary of the inviscid flow is shifted is called the bound- 
ary layer displacement thickness, and is denoted by 8 . The displacement thick- 
ness can be either positive or negative, depending on the boundary layer velocity 
and temperature profiles. In typical NATA solutions, 8 * is negative in the up- 
stream and throat regions, and becomes positive at some point downstream of the 
throat. In channel flow problems, there are two displacement thicknesses, 8 i 
and 82 * , one for each pair of channel walls. These displacement thicknesses are 
generally unequal because of the different amounts of streamline convergence or 
divergence in the boundary layers on the two pairs of walls. 


In many parts of the code, it is necessary to calculate the geometric area 
ratio Ag from the effective area ratio A e , or conversely. The relation between 
the two area ratios, upon which both of these subroutines are based, depends upon 
the type of nozzle geometry. For a two-dimensional nozzle, the effective flow 
area per unit length in the 2 -direction is 

Ag = 2 [y (x) - 8 *] (124) 


At the sonic point, the effective area is 

A ' = 2 [y* - 5*] 

* 


(125) 


where the subscript * denotes sonic conditions. Hence, the effective area ratio 
is equal to 


A 


e 




A' 


* 


y s * 
y -8* y* y* 



y* 


(126) 
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Because the sonic point is near the geometric throat, and because dy/dx = 0 at the 
throat, the profile ordinate y* at the sonic point is approximately equal to the 
ordinate yo at the throat. In NATA, the difference between these two quantities 
is neglected. Hence, from (118) and (126), 



For an axisymmetric nozzle, the effective flow area is 

A^ = ff (y-5*) 2 (128) 

and at the sonic point 


A ' = n (y* - S*) 2 

Hence, with the approximation y * — yq 


(129) 



Equations (116) and (130) give 



(130) 


(131) 


In the case of a rectangular channel, the effective flow area is 

Ag = 4(y - Si) (z - S* 2 ) 032) 

while at the sonic point 

A' = 4.(y*-Si*)(z»-4) (133) 

* 

Hence, with the approximations y* ^yQ > z * ~ z o » 
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( 135 ) 


Combination of (120) and (134) gives 


A 


g 


1 8 Ma yS 2 + z8 1 

~ To ~*oj e+ yo *0 



70 2 o 


4 * 

This is not an explicit solution for A g as a function of A e and s i , S 2 , be- 
cause the right hand side still contains the profile ordinates, y andz . These 
are known as functions of * , but * is not known until A g has been determined. 

It is therefore necessary to carry out the solution for A g using an iterative 
method, in the case of channel geometry. 
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5. LAMINAR BOUNDARY LAYER 


The NATA code contains an optional calculation of the buildup of a laminar 
boundary layer on the nozzle or channel wall. The results of this calculation 
are used in two distinct ways: 

(1) The boundary layer displacement thickness S is used to relate the ef- 
fective inviscid area ratio A e to the geometric area ratio A g , as ex- 
plained in Section 4. In this way, the inviscid flow solution is 
coupled with the boundary layer solution. 

(2) The heat flux and shear stress at the wall are printed as code outputs. 
These data can be used to test the overall flow model by comparison with 
experimental measurements. In the case of a rectangular channel con- 
taining a test panel, these outputs provide predictions of the test con- 
ditions on the panel. 


5.1 Boundary Layer Transition Criterion 

The laminar boundary layer calculation in NATA is applicable, of course, only 
if the boundary layer in the case being investigated is actually laminar and not 
turbulent. Schlichting (ref. 21, pp. 457-463) has reviewed the available data on 
boundary-layer transition in incompressible flows. A stability analysis shows 
that a laminar boundary layer is unconditionally stable as long as the Reynolds 
number based on momentum thickness, 

N Re<? = Pe u e (136) 

is less than 162. In (136) , p e , , and n e are the free-stream density, velocity, 

and viscosity, respectively, and 6 is the momentum thickness. 



(137) 


In (137), y is a coordinate locally normal to the surface. For values of N Re $ 
higher than the critical value of 162, the boundary layer is unstable, but the 
occurrence of transition depends upon the level of turbulence in the free stream. 
Transition occurs near N Re g = 162 only for extremely high turbulence levels. For 
exceptional] y smooth free-stream flows transition does not occur until an N Re ^ of 
about 90C is reached. 


the transition Reynolds number also varies with other flow parameters. It 
increases with increasing Mach number and with’ increasing favorable pressure 
gradient, and decreases with mass injection at the wall. In a recent study (ref. 
22), extensive data on boundary layer transition have been compiled and correlated 
with free-stream Mach number M e and N Re $. The data can be fitted by the relation 


N Re0, T “ 


0.224 M 

200 e 


e,T 


(138) 
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where the subscript T refers to conditions at transition. The data contain con- 
siderable scatter, so that they do not determine the numerical coefficients in 
(138) with high precision. The choice of the value 200 for the critical Reynolds 
number at zero Mach number, in place of 162 or other possible values, takes ac- 
count of the fairly high level of turbulence present in arc-heated air streams, 
of the strong favorable pressure gradient accelerating the flow in a converging- 
diverging wind tunnel nozzle, and of the absence of mass injection through the 
nozzle wall. At Mach numbers below 6, equation (138) gives a lower transition 
Reynolds number than the correlation recommended in reference 22 for flight 
conditions . 

The Reynolds number Nj^ e g (136) and the transition Reynolds number q- (138) 
are computed and printed out by NATA. In most cases lying within the operating 
envelopes of existing NASA/JSC arc-heated wind tunnels , N Re $ remains smaller than 
^Re0 , t throughout the solution.* Thus, the assumption that the boundary layer is 
laminar is probably valid in most cases of current interest to NASA/JSC. However, 
a proposed upgraded pumping capability allowing mass flows up to 1.5 lb/sec would 
make accessible a range of operating conditions in which the boundary layer would 
often become turbulent. 


5 . 2 Basic Equations and Transformations 

The basic non-equilibrium inviscid flow solution carried out by NATA requires, 
typically, a few minutes of computer time. The equilibrium and frozen inviscid 
solutions are much faster. The choice of a method for performing the laminar 
boundary layer calculation was guided, in part, by the requirement that solutions 
including the boundary layer should not consume a great deal more computer time 
than the basic inviscid solutions. This criterion ruled out the use of exact 
solutions of the partial differential equations for the laminar boundary layer, 
and narrowed the choice to one of several available approximate methods. The 
technique chosen is an integral method developed by Cohen and Reshotko (ref. 23). 
This technique is an extension of Thwaites' correlation method for incompressible 
boundary layers (ref. 24) to the compressible case with heat transfer. It treats 
the effects of arbitrary variable pressure gradients. According to Hayes and 
Probstein (ref. 25, pp. 318-320) the Cohen-Reshotko integral method is probably 
superior in accuracy to the method of local similarity. 

The Cohen-Reshotko method uses curvefits to certain boundary layer properties 
based on similar solutions. In their original report (ref. 23), Cohen and 
Reshotko based these curvefits on similar solutions for Prandtl number unity and 
viscosity proportional to the absolute temperature. The corresponding curvefits 
used in NATA are based on similar solutions carried out by Dewey and Gross (ref. 
26), and include the dependence upon Prandtl number, the viscosity-temperature 
index <u appearing in y. « 1*°, and the hypersonic parameter a = u e 2 / 2hg. Thus, there 
is reason to hope that the boundary layer calculations in NATA may be somewhat 
more accurate than those given by the Cohen-Reshotko method in it's original pub- 
lished form. 


*An exception is the boundary layer on the non-expand|ng face of a rectangular channel, which often becomes turbulent a short 
distance beyond the throat. 
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The relations used in the Cohen-Reshotko method will be derived starting 
from the equations for an axisymmetric boundary layer in a compressible fluid. 
(See ref. 27, pp. 384-385.) 




dp 

dx 



(139a) 


d d 

— (pur) + — (p vr) 
dx dy 


(139b) 


dh 


dh 


pu— + pv — - 
dx ay 


dp 

dx 


d/p dh\ /d u \ 2 

dy \N Pr <?yj \dy j 


(139c) 


In these equations, p denotes the density, p the pressure, P the viscosity, b the 
static enthalpy, and Np r the Prandtl number. The independent variables, x and y , 
are a pair of locally orthogonal coordinates, with x parallel to the surface along 
the streamwise direction and y lying along the normal to the surface. The local 
radius of the surface from the axis of symmetry is denoted by r(x). 


The equations (139) for an axisymmetric boundary layer can be converted into 
those for a two-dimensional boundary layer using the Mangier transformation (refs. 
21 and 27): 


// 

*0 


r 2 (x) dx 


r(x) 


yu dr 


v = — v + 

r \ r dx 


(140a) 


(140b) 

(140c) 

(140d) 


where L is an arbitrary non-zero scale length. The equations for a two-dimensional 
boundary layer are 


pu 


du 

dx 


_ du dp d / du \ 

+ pv + (p 
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(141a) 

(141b) 

(141c) 
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Cohen and Reshotko (ref. 23) base their analysis on equations (141). They 
first convert these equations into the equations for an incompressible boundary 
layer using the Stewartson-Illingworth transformation (refs. 21, 28, 29): 


X = 



a e (x ') p e (x) 

dx ' 

**0 PQ 


(142a) 


Y = 


a e (x > f p(x, y') 
I dy 

a 0 J P0 
J o 


(142b) 


Here, 


A (SO 


Pq/Tq 


(143) 


where T represents the absolute temperature and the subscripts w and 0 denote con- 
ditions at the surface (wall) and free-stream stagnation conditions, respectively. 
The symbol a represents the sound speed, and subscript e indicates local free- 
stream (external) conditions. The transformation of the velocity components is 
found by defining a stream function 'A such that 

d\j/ p u 

dy p 0 

dijj p v 

dx pq 

Then 


(144a) 

(144b) 


U 


dtfr 

Jy 



v - OY /dx)- (ag/a^u 
dX HpQ/p) (ag/ag) (Pj/po) 

To obtain the desired results it is necessary to assume that 


(145a) 


(145b) 



( 146 ) 


- 54 - 



which is equivalent to 



( 147 ) 


in which A is given by (143). It is .also necessary to assume that the inviscid 
flow obeys the equations for a perfect gas, in particular 



(148a) 


(148b) 


The transformed momentum and continuity equations for the boundary layer are then 


U 


<3U 

lx 


+ 




+ 


(1 + S) u e 


du e 

dX 


au 

ax 



(149a) 

(149b) 


Cohen and Reshotko also give the transformed energy equation, but this is not re- 
quired in the subsequent analysis. In equation (149a) 



(150) 


where h s is the local stagnation enthalpy and ho is the stagnation enthalpy in the 
free stream, and 

"o = «/po 

The transformed boundary conditions are 


U (X, 0) = 0 (152a) 

V(X, 0) = 0 (152b) 

S (X, 0) = S W (X) (152c) 

lim S = 0 (152d) 

Y CO 

lim U = U e (X) (152e) 

Y -> oo 
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5.3 Momentum Integral 


The incompressible momentum integral equation is now obtained by multiplying 
the continuity equation (149b) by (U e -U) , and subtracting the resulting relation 
from the momentum equation (149a). This gives 


i [U(U, - U)1 + 1 tvro. - TO] 


du, 


e 


+ SUe lx + v ° 


dY z 


(153) 


This equation is now integrated with respect to Y from Y ■ 0 to A, where A is a 
fixed distance encompassing the entire boundary layer. The result is 




(Si* + 6) 



where 



(154) 


(155a) 


(155b) 


(155c) 


are the incompressible momentum, displacement, and enthalpy thicknesses for the 
boundary layer, respectively. Equation (154) is equivalent to 




(2 ^ + 8 * + g) 


dU„ 


dX 



(156) 


- 56 - 



Cohen and Reshotko define the following non-dimensional boundary layer 
parameters : 


Shear parameter 

u e 


(157) 


*i 2 

du e 


Correlation parameter 

n = — ■ — — 

v 0 

dX 

(158) 


In terms of these parameters, the equation (156) can be expressed (after multipli- 
cation with 2 0j U e /vq ) in the form 


— U. 


dX 


(-3u73f) " 2 (H »c + 3 * *1 


(159) 


where 


H:, 


8 * + € 


is termed the incompressible form factor. 


(160) 


5.4 Correlation Method 


Cohen and Reshotko now assume that the non-dimensional boundary layer param- 
eters l and H; nc , and also the Reynolds analogy factor (to be defined later) , can 
be approximated as functions of the correlation parameter n and the value S w of S 
at the wall, irrespective of the previous history of the boundary layer. This is 
their generalization of Thwaites’ correlation concept. On this assumption, the 
right-hand side of equation (159) is a function of n and S w : 

N(n, S w ) = 2 [n (H inc + 2) + £] (161) 

so that (159) becomes 


CL 



N(n, S w ) 


(162) 


Figure 1 shows the form of the function N(n, S w ) for favorable pressure gradients 
(n<0 ) , based on similar boundary layer solutions for Prandtl number unity and 
the viscosity law (146). The points in this figure represent results tabulated by 
Dewey and Gross (ref. 26). The parameters n ,£ , H inc , and N appearing in Cohen 
and Reshotko' s analysis are related to the quantities f"(0) , 1^ , I 2 » and P tabu- 
lated by Dewey and Gross as follows: 


n = ~ P I 2 2 


(163a) 
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Z = I 2 f "(0) 


H- 1 

me t 

h 


(163b) 

(163c) 


N = 2I 2 [f"(0) - j8(I 1 + I 2 )]' 


(163d) 


The family of straight lines in figure 1 represents a simple curvefit to the Dewey 
and Gross Data: 


N (n, S w ) = A + B (S w ) n 


where 


B(S W ) . 2.14 + 1.28 (1 + S W ) + 0.93 (1 + S W ) 4 


(164) 


(165a) 


and 


A = 0.38 


(165b) 


The linear dependence of (164) upon n is required for further progress in the 
analysis. The figure shows that the linear fit is reasonably accurate. At n = 0, 
the correct value of N (n t S w ) is 0.441 for all S w , whereas the analytical fit 
(164), (165) gives a value of 0.38. The coefficients (165) have been chosen to 
optimize the accuracy for n <_ 0.1, i.e., for rather strongly favorable pressure 
gradients, because this is the region of most frequent use and greatest importance 
in NATA calculations. 

Substitution of the curvefit (164) into the differential equation (159) gives 


-a 


d 

e dX \ dU VdX 


= A + B(S W ) n 


(166) 


For constant wall temperature (S w = const.), equation (166) has the solution 


n = - AU. 


dX 


7 u 

\ 


B- 1 


dX 


(167) 


where Xq is the value of X at the point where the boundary layer starts. Initial 
conditions on the boundary layer will be discussed below. 

The external flowfield U e (X) is assumed known. Thus, n can be calculated at 
each point X by a running quadrature, using (167). Then the incompressible 
momentum thickness can be obtained from (158) : 






dUg/dX 


( 168 ) 
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Because the shear parameter £ is assumed, to be given as a function of n and S w , 
and because S w is known, the shear can then be calculated from equation (157). 
Similarly, the displacement thickness can be evaluated from Hi nc (n,S w ), and the 
heat flux to the wall from a Reynolds analogy factor. Thus, equation (167), to- 
gether with Cohen and Reshotko's correlation assumption, provides an explicit 
description of the boundary layer. 

5.5 Momentum Thickness 

Substitution of the Stewartson- Illingworth transformation formulas (142a) 
and (145a) gives the solution (167) for n in terms of the variables for a com- 
pressible two-dimensional boundary layer: 


n 



(169) 


Then substitution of the Mangier transformation formulas (140a) and (140c) gives 
the expression for n in the case of an axisymmetric compressible boundary layer: 


n 


A M e -B 

2 3e If. 

a 0 Po 



(170) 


To introduce a more compact notation, let 

<D . r 2 > Mj 3 - 1 — — (171) 

*0 PO 

I = J $ (172) 

0 


a' = Al/L 4> (173) 

where L is a characteristic length and the index j is 0 for a two-dimensional 
boundary layer and 1 for an axisymmetric layer. Then 


n 


Ln' 



dM e 

dF 


(174) 
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t 


where £ is the streamwise coordinate along the surface (equal to x - xq for the 
axisymmetric case and to s - ^ for the two-dimensional case) . Both n ' and n are 
non-dimensional . 

The momentum thickness for a compressible two-dimensional boundary layer is 
defined by 


f 8 

/ _£L (l-~ 
J Pe u e \ u e 


dy 


(175) 


Equations (142b), (145a), (155a), and (175) show that 6 is related to the incom- 
pressible momentum thickness by 


Pq a 0 
Pe a e 


(176) 


Application of the Stewartson-Illingworth formulas (142a) and (145a) to equation 
(168) , followed by use of (176) , gives the following explicit formula for the 
momentum thickness of a two-dimensional compressible boundary layer: 


Q = L 


n M e 

1 X 

PO 

PO P e 

dM e / dx 

N ReL L 

15 

n> 

Pe PO 


(177) 


where 


N ReL = 


Pe u e L 

Pe 


(178) 


is the Reynolds number based on the characteristic scale length L . The value of 
0 is independent of L . Because 4 is defined by (143) , and since the pressure 
gradient through the thickness of the boundary layer is negligible (P w = P e ) , for 
a perfect gas 


PO PO P e Pw Pw 

Pe Pe PO Pe Pe 


Hence, (177) becomes 


j nM^ i 

y ~ dMg/d x ’ L Re L 


Pw Pw 
Pe Pe 


The momentum thickness for an axisymmetric boundary layer is 

.8 


e 


f _£L_ (i _ JL\ 

J Pe u e \ u e / 


dy 


(179) 


(180) 


(181) 


- 61 - 



From the Mangier transformation (140) , 

e = — 0 (182) 

r 

and 


d M e L 2 dM e 

dx" r 2 dx 

Substitution of these relations into (180) gives 


e = l 


j ° 

dMg/dx 


Pw Pw 


u e/ d* LN ReL Pe Pe 


(183) 


(184) 


which is of the same form as (180). Thus, based on this formulation, the differ- 
ence in momentum thickness between axisymmetric and two-dimensional boundary- 
layers with the same external flow M e ( f ) arises entirely from the difference in 
n resulting from use of equation (169) in the 2D case and (170) in the axisym- 
metric case. In both cases, the momentum thickness can be expressed, with the 
aid of (174) , in the form 


Q = L 



Pvf Pw 
Pe Pe 


(185) 


The above formulation assumes, in the axisymmetric case, that the boundary 
layer thickness is everywhere negligible compared with the distance r of the wall 
from the axis of symmetry. This assumption is not necessarily valid far down- 
stream in an axisymmetric nozzle, where the local Reynolds number N ReL can become 
quite low owing to the rarefaction of the gas. Accordingly, for axisymmetric 
nozzles, NATA applies an approximate transverse curvature correction to the momen- 
tum and displacement thicknesses. This correction is based upon the following 
reasoning (reference 30). Figure 2 shows a portion of an axisymmetric nozzle, in 
cross section, near a location where the nozzle radius is r and the tangent to 
the nozzle surface is inclined at an angle b to the nozzle axis. The actual momen- 
tum thickness 6' , including the effects of transverse curvature, is shown. The 
flow area in the boundary layer, based on the thickness 6' , is the area of the 
conical frustum whose cross section is represented by the line segment AB in the 
figure. This area is given by 


A = n (2 r — 6 ' cos b) 6 ' 


(186) 


It is assumed that this flow area is equal to the flow area 2ncd, where 6 is the 
momentum thickness (184) or (185) based on the thin-layer analysis. This assump- 
tion gives the relation 



( 187 ) 


This correction increases the calculated momentum thickness for axisymmetric 
nozzles. 
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Figure 2 GEOMETRY OF TRANSVERSE CURVATURE CORRECTION 
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5.6 Displacement Thickness 


Cohen and Reshotko (ref. 23) give the following relation between the displace- 
ment and momentum thicknesses: 


8 * 

-J - (Hinc + » 



(188) 


where Hj nc is defined by (160) . This relation can also be derived from relations 
given by Dewey and Gross (ref. 26) for similar boundary layers with arbitrary 
Prandtl number, free-stream Mach number, and viscosity-temperature index to , for 
zero sweep angle. NATA calculates 8 from (188), using an analytical curvefit to 
Hinc based on the Dewey-Gross data. For axisymmetric boundary layers, NATA uses 
the momentum thickness 6' corrected for the transverse curvature effect, equation 
(187), in (188). 

The curvefit to the incompressible form factor is 


H ; nc + 1 = {0.25 +0.75 e 4 - 6n [l + (3- 15-6n)(l+S w )]t • 


td + S w ) 


-(M&iNp, 


-0.47 a In Npj] 


(189) 


where Np r represents the Prandtl number and or is the hypersonic parameter 



(190) 


This formula contains two main factors. The first, enclosed in braces, gives the 
dependence of Hi nc upon n and S w for Prandtl number unity and = 1. The second 
factor, in brackets, gives the dependence on Np t and a . The dependence of Hj nc 
upon the viscosity-temperature index o> , defined by 




(191) 


is weak and is neglected in (189). Figures 3 to 5 compare the curvefit (189) with 
data on H^c calculated from Dewey and Gross’s tabulation of similar solutions. 

The abscissa in these plots is 1+ Hj nc as calculated from (189) , while the ordinate 
is 1 + Hi nc based on the Dewey and Gross data. The plots include data for Np r = 0.5, 
0.7, and 1.0, for <o = 0.5, 0.7, and 1, for o = 0, 0.5, and 1, for (1 + Sw) = T w /T 0 
from 0 to 1, and for values of jS (the Falkner-Skan parameter) from 0 to 5. The 
root-mean-square percentage error in (1+Hj nc ) based on the curvefit, for the 405 
data points considered in the fitting analysis, is 5 percent. The largest in- 
dividual error is 17 percent. 

In the limiting case of an infinite favorable pressure gradient (jS-»°o) , Dewey 
and Gross's integral Ii is zero (reference 26). Thus, from (163c), Hjjjj.-^-l. The 
curvefit (189) does not reproduce this limit exactly. For ^-»°°,n->-oo according to 
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Figure 3 COMPARISON OF INCOMPRESSIBLE FORM FACTOR CURVEFIT WITH 
DEWEY-GROSS DATA (N Pr = 0.5) 
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Figure 4 COMPARISON OF INCOMPRESSIBLE FORM FACTOR CURVEFIT WITH 
DEWEY-GROSS DATA (Np r = 0.7) 
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Hj nc + 1 (CURVEFIT) 

Figure 5 COMPARISON OF INCOMPRESSIBLE FORM FACTOR CURVEFIT WITH 
DEWEY-GROSS DATA (N Pr = 1.0) 
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(163a); thus, the term containing the factor e^ n in (189) has a limit of zero, 
and the limit of ( 1 + Hj nc ) is 

lim [1 + H inc (fit)] = 0.25 [(1 + S w f°- 4 ^ Npr - 0.47 <7 £ n N Pr ] ( 192 ) 

-> oo 

For example, if Np r = 0.7, (1 + S w ) = 0.1, and a = 1, this limit has the value 0.22. 
The correct value for j8-»~ would be zero. However, the Dewey-Gross data show 
that the approach to this limit is very gradual. For the pressure gradients en- 
countered in NATA runs, (189) is of acceptable accuracy. 


5.7 Shear Stress 

The shear stress at the wall can be calculated from the definition and 
value of the shear parameter £ . From equations (157), (142b), (145a), (140b), 
and (182) , 


Pw f*w l “e 
Pe 6 


(193) 


for either an axisymmetric or a two-dimensional boundary layer. 


Values of £ can be obtained from Dewey and Gross's tabulated similar solu- 
tions using equation (163b). Values of the shear parameter thus obtained have 
been curvefitted by the following expression: 


£ = 


0.2205 + 


1 + 4.5(1+S W ) 0 - 9 


2.7 


(-n) 


5.06 

(_n)°-224 


<£o + o.45 (1 - co) 




where 


<b = exp | - 0.427 [(l + S w ) -(1 -' a) -l] } 


^1 


0.76 (1-cu) 1 - 2 (1 - VTT’a) 


l + a(-n) b 


a = 15.5 e 


- 29(1+S w> + 0.67e 


4.37(1 + S„) 


(194a) 


(194b) 

(194c) 

(194d) 


b = 0.7 + 0.47 (1 + S W ) 


(194e) 


The first factor in (194a) gives £ as a function of n and S w for ® *• 1. The 
second factor gives the dependence on <a and o . The dependence on Prandtl number 
is weak and is neglected. Thus, equations (194) represent £ as varying with the 
hypersonic parameter a even for Np r =1. In theory, the boundary layer properties 
should be independent of or fbr Np r = l. The Dewey-Gross data indicate little 
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dependence of t on Np r in the range from Np r = 0.5 to 0.7, and thus provide no 
basis for fitting a Prandtl number dependence which would cause the dependence 
on o to vanish for Np r = l. These formulas should, in any event, be satisfactory 
for air, in which the Prandtl number is about 0.7. 

In typical applications of NAT A, the enthalpy ratio parameter ( 1 + S w ) is 
small of order 10“2. The smallest value of (1+^) in the Dewey-Gross data, for 
<u <1, is 0.05. Thus, application of these data in the code normally requires an 
extrapolation to (1+S W ) values beyond the range investigated by Dewey and Gross. 
Unfortunately, the data for (1+^) = 0.05 and 0.1 indicate that £ has a strong 
dependence on (1+S W ) for small (1 + S w ) and w <1. As a result, the extrapolation 
to still smaller values of (1 + ^) is a potential source of sizeable errors. The 
rather complex curvef it (194) was developed in an effort to minimize such extrapo- 
lation errors by fitting the Dewey-Gross data for small (1+S W ) as closely as 
possible. 

Figures 6 through 8 compare the curvefit (194) with l values obtained from 
Dewey and Gross's tabulated similar solutions. The 405 cases used in the compari- 
son include Prandtl numbers of 0.5, 0.7, and 1, <u values of 0.5, 0.7 and 1, values 
of the hypersonic parameter a of 0, 0.5, and 1, several temperature ratios 1 + S w 
in the range from 0 to 1, and values of the Falkner-Skan parameter /3 ranging from 
0 to 5. The root-mean-square percentage error is 2 percent. The largest individual 
percentage error is 17 percent. 


5. 8 Heat Flux 

Because the energy equation is not used in the Cohen-Reshotko correlation 
method, the heat flux is calculated from the shear stress r w using a Reynolds 
analogy factor based upon similar solutions. The formula is 

r w R A u e N Prw 

where h r denotes the recovery enthalpy and Ra the Reynolds analogy factor, 
expression for recovery enthalpy recommended in the literature (ref. 25, p. 
is 

h r = h e + VNp r (h 0 - h e ) 

However, the data tabulated by Dewey and Gross are fitted more accurately by 

h r = h e + N Pr 0 - 56 (h 0 - h e ) (197) 


(195) 

The 

296) 

(196) 


which is the formula used in NATA. The Reynolds analogy factor is given, in terms 
of quantities tabulated by Dewey and Gross, by 


*A = — 
A d'(0 ) 


£ aw t w 



(198) 
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Figure 6 COMPARISON OF SHEAR PARAMETER CURVEFIT WITH DEWEY-GROSS DATA 

(N Pr = 0.5) 
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Figure 7 COMPARISON OF SHEAR PARAMETER CURVEFIT WITH DEWEY-GROSS DATA 

(N Pr = 0.7) ' 
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Figure 8 COMPARISON OF SHEAR PARAMETER CURVEFIT WITH DEWEY-GROSS DATA 

(N Pr = 1.0) 
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( 199 ) 


Here, t aw = h r /ho . A curvefit to t aw , based upon equation (197), is 

*aw = 1 - a ( 2 ~ N Pr 0 - 56 ) 

The Reynolds analogy factor based upon Dewey and Gross's tables has been curve- 
fitted in the following form: 


R a = N, 


-0.35 


Pr 



n 

a R ~ b R n 


(200a) 


where 

a R = (l + S^" 0 - 86 " 

+ (1 -N Pr ) 2 [-0.185(1 + 2 S w )+ 0.29 (l+S w )°- 38 a/o>] | (200b) 

b R = 0.94 - 1.31 (1 + S W )°- 486 


6.2 + 2.35 (l + SJ- 0 - 86 


- in Np r (0.62 (1 + S w ) + [3.5 (1 - a) - 0.49 - 3,61 (1 + S w )] a | 


(200c) 


As in the case of the shear parameter, a relatively complex analytical form is 
used in an attempt to minimize the errors resulting from extrapolation to values 
of (1+S W ) below the range covered by the Dewey-Gross data. 

Figures 9 through 11 compare the curvefit (200) with data from Dewey and 
Gross. The ordinate in these figures is B'(0) /f"(0) , computed from values in 
their tabulation of similar solutions. The abscissa is the quantity 
[(t aw - **)/(! ~ c w^/ r A (00 > calculated from equations (199) and (200). Equation 

(198) shows that the ordinate and abscissa would be equal if the curvefit were 
exact. This manner of plotting the comparison tests the curvefit (197) for the 
recovery enthalpy together with the fit (200) for the Reynolds analogy factor. 
Figures 9 through 11 do not include the Dewey-Gross points for which t w = 1 or 
t w “ W • Such points were excluded because it was difficult to fit both the 
high-wall-temperature and low-wall-teraperature heat transfer data accurately using 
a single analytical formula. The fit to Ra thus applies only to cases in which 
the wall temperature is lower than the adiabatic wall temperature, which is always 
the case in typical applications of NATA. The plots in figures 9 through 11 con- 
tain a total of 321 points. The root-mean-square percentage error for these 
points is 9 percent. The largest individual percentage error is 127 percent, in 
a case with t w = 0.6 and t a w = 0.67. The error arises partly from inaccuracy in 
the curvefit to the adiabatic wall temperature. Another similar point has an 
error of 79 percent. If these two points were excluded, the maximum percentage 
error would be 24 percent and the root-mean-square percentage error would be 
3 percent. 
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Figure 9 COMPARISON OF REYNOLDS ANALOGY FACTOR CURVEFIT WITH 
DEWEY-GROSS DATA (Np r = 0.5) 
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Figure 10 COMPARISON OF REYNOLDS ANALOGY FACTOR CURVEFITWITH 
DEWEY-GROSS' DATA (N Pr = 0.7) 
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Figure 1 1 COMPARISON OF REYNOLDS ANALOGY FACTOR CURVEFIT WITH 
DEWEY-GROSS DATA (N Pr = 1.0) 
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5.9 Dependence of the Momentum Parameter N Upon Np r and co 


Cohen and Reshotko assumed the parameter N, appearing in the form (162) of 
the boundary layer momentum equation, to be a function only of n and S w . When 
N is evaluated from the Dewey-Gross data using (163d) it is found to vary with 
Np r , co , and a as well as with n and S w . Cohen and Reshotko neglected the dependence 
upon these additional parameters because they determined N from similar solutions 
for Np f = co = 1. For these conditions, the boundary layer properties are inde- 
pendent of a . 

Because the solutions tabulated by Dewey and Gross include the dependence 
upon Np r , co , and a , it appears desirable to determine whether this information 
can be used to improve the accuracy and range of applicability of the Cohen- 
Reshotko correlation method. The solution (167) of equation (166) is valid only 
if the coefficients A and B in the linear curvefit (164) are constant for each 
problem. Since a varies within each problem, from zero in the upstream reservoir 
to nearly 1 at the downstream nozzle exit, including the dependence on o would 
lead to inadmissible variation of A and B. Thus, the dependence on o must be 
neglected. 

Inclusion of the dependence of N upon the Prandtl number leads to no incon- 
sistency or difficulty so long as Np r is constant for each problem. However, in- 
clusion of the dependence upon <o involves a logical inconsistency, because the 
derivation of equation (162) involves the assumption of the viscosity law (146) , 
which is equivalent to (191) with co = 1. Nevertheless, a dependence of A and B 
upon co is compatible with the solution (167), so long as co is constant for each 
problem and the inclusion of such a dependence can be justified by comparison with 
Dewey and Gross's similar solutions (see below). 

The variation of N with N Pr and co can be fitted approximately by retaining 
the analytical form (164) and the expression (165b) for B, but replacing (165a) by 

-6.67 (1 + S ) 

A = 0.38 — 0.76 Np f (1 — <u) e w (201) 


Figures 12 through 14 compare the curvefit given by (164) , (165a) and (201) with 
the Dewey-Gross data. The straight line in the figure is the locus of agreement 
between the fit and the data. The curvefit represents the data reasonably well 
for N values less than about 0.2, but is systematically too low at higher N . 

This is the same discrepancy seen near the right-hand edge of figure 1, where the 
linear curvefit N = A + Bn lies beneath the upward-curving relationship N ( n , S w ) 
given by the Dewey-Gross data. The region of N values where the curvefit is most 
accurate corresponds to values of the correlation number n less than about -0.05, 
i.e., to the range of strongly favorable pressure gradients which is important 
in NATA runs for axisymmetric nozzles. For channels, smaller values of n of 
order -0.01 are encountered on the expanding face. 


The validity of treating the parameter A in equation (166) as a function of 
co and Np r , in accordance with (201), can be tested by applying the Cohen-Reshotko 
integral method to some of the similar solutions tabulated by Dewey and Gross. 
These similar solutions assume that the Falkner-Skan parameter /3 is constant. By 
definition 


2 £ du e T 0 

u e d £ T e 


( 202 ) 
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Figure 12 COMPARISON OF MOMENTUM PARAMETER CURVEFITWITH 
DEWEY-GROSS DATA (N Pr = 0.5) 


- 78 - 






DEWEY-GROSS DATA (Np r = 0.7) 
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Figure 14 COMPARISON OF MOMENTUM PARAMETER CURVEFIT WITH 
DEWEY-GROSS DATA (Np r = 1.0) 
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(203) 


■/ 


PwPw u e 


r 2 J dx 


For p w ^ = constant and j =0, (203) becomes 


f = 



dx 


(204) 


To simplify the analysis, cases with a « 1 are considered. Then (202) gives, 
upon integration with To/T e = 1, 


f = K Ug 2 // 3 

(205) 

With the aid of (204) and 
ential equation with x as 

(205) , equation (202) can be converted into a dif fer- 
tile independent variable: 

2K u e (2//3)- 2 du e 

r) t r 


Pw Pw P 

(206) 

For a « 1, 


a e 

— r 1 
*0 

(207a) 

Pe 

- 1 

Po 

(207b) 

M e - Ug/ag 

(207c) 


Thus, for the cases under consideration, the quadrature in equation (169) can be 
performed analytically with the aid of (206) and (207) . The result is 


A 

B — 2 + (2//3) 


(208) 


With equation (164) , this can be written in the form 


N (fit) 


(209) 
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Figure 15 compares values of n calculated from (209) with values obtained 
directly from (163a) , for a number of the similar solutions tabulated by Dewey 
and Gross. The filled circles represent data for &> = 1, Np r = l. The other 
symbols refer to cases with co A 1. Cases with >3=1, for which the denominator in 
(209) is zero and N(fit ) is small, are omitted. The line is the locus of agree- 
ment. The figure shows that the n values calculated from (209) agree reasonably 
well with the correct values based on (163a) . The agreement is nearly as good for 
co U 1 as for (o = l. If the dependence of A upon co and Np r , as given by (201), 
had not been included in the calculation, the agreement would have been poorer 
for the cases with co 4 1, especially for smaller-magnitude values of n . Thus, the 
use of (201) to determine the coefficient A in equation (166) appears justified. 


The data in figure 15 are all for o = 0. If data for other values of o had 
been included, the scatter would have been greater, especially for small |n | . 


It is also of interest to compare the predictions of momentum thickness based 
on the Cohen-Reshotko integral method with the exact results for the similar solu- 
tions. Substitution of equations (204) to (208) into the formula (177) for the 
two-dimensional momentum thickness gives, with use of (179), 
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( 210 ) 


The momentum thickness for similar solutions, based on relations given by Dewey 
and Gross, is 
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( 211 ) 


where £ has been expressed using (205) . From (210) and (211) , 
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According to equations (163a) and (208) , the quantity under the radical is simply 
the ratio of the correlation number n , as evaluated using the Cohen-Reshotko 
method with the curvefit (164), to the correct value of n for the similar solution. 
Thus, to the extent that (208) agrees with (163a), the Cohen-Reshotko integral 
method gives the correct momentum thickness when applied to similar boundary layer 
problems, at least in the limit a « 1 that is treated by the preceding analysis. 

The remaining boundary layer properties ,5* , r w , and q w , are all calculated 
from the momentum thickness using curvefits to relations based on similar solutions. 
Thus, their accuracy should be comparable with that of the momentum thickness. 


5.10 Initial Condition 

The initial condition on the solution (167) of the Cohen-Reshotko differential 
equation (166) is embodied in the lower limit of integration, Xq . If U e is 
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Figure 15 COMPARISON OF THE CORRELATION NUMBER n (CALCULATED USING THE 
COHEN-RESHOTKO APPROXIMATION) WITH DATA FROM SIMILAR SOLUTIONS 
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non-zero at X = Xo , then n = 0 at that point. In that case, the incompressible 
momentum thickness is also zero at X Q , If X = Xq is a stagnation point 
(Ug = 0) then n(Xo) is not zero and the boundary layer starts with a finite 
thickness . 

It is difficult to specify a fully satisfactory initial condition for the 
boundary layer in an internal flow system such as an arc-heated wind tunnel. 

Under some circumstances the boundary layer might extend upstream along the wall 
confining the flow, through the plenum chamber (if any) and the arc heater. The 
conditions in these regions upstream of the nozzle are poorly known and generally 
chaotic because of the presence of the arc column, swirl caused by tangential 
gas injection and, in some cases, free-stream turbulence. In general, the cross 
sectional area of the flow is finite everywhere , and any stagnation point which 
may exist in a "corner" of the wall profile is likely to receive an influx of 
boundary layer gas from further upstream. 

In NATA, no attempt is made to model the detailed geometry of the plenum and 
arc chamber. The boundary layer is assumed to start at a somewhat arbitrary point 
near the juncture between the nozzle and the plenum. The uncertainty in x 0 due to 
the arbitrariness in the choice of this starting point affects the value of the 
integral I (172). The resulting uncertainty in I, in turn, affects n' (173), 
n (174), and thus 0 (185), 8* (188), (193), and q w (195). The percentage errors 

in these boundary layer properties, due to an error in xq , decrease downstream as 
the value of I increases. 

Some insight into the dependence of these errors upon x can be obtained by 
examining the behavior of the integral I in the limits of low and high Mach number. 
At the assumed boundary layer starting point xQ , the Mach number is much less than 
unity. Thus, a e ~ Pe/PO ~ 1 , and I is approximately 



In the two-dimensional case ( j = 0) , the contributions to I from points near £ = 0 
are relatively small because u e is low in that region and, according to (165a), 

B -1 is a little greater than 1 . Thus , for two-dimensional boundary layers , the 
results in the throat region and downstream should be rather insensitive to the 
value of xq chosen, so long as it defines a point well upstream of the throat 
where the Mach number is low. 

In the axisymmetric case (j = 1), however, the continuity equation pu e A'= 
constant shows that r 2 i u e B ~ is roughly constant in the subsonic part of the 
flow, because r 2 ) is proportional to the cross sectional flow area A', 3~lis 
approximately 1 , and p is almost constant. Thus, the results for an axisymmetric 
boundary layer are more sensitive to the choice of xq than those in the two- 
dimensional case. 

For flow over a flat. plate at zero angle of attack, the integrand $ (171) of 
I is constant, so that I <* (x-xq ) . In that case, the relative error in I due to 
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an error in xg decreases as (x-xg) For flow of a perfect gas in an axisymmetric 
nozzle, the integral can be expressed in the form 


I = Rg 2 



0.5(y+l)/(y- 1) 




M, 


B — 3 




(214) 


In the limit of high Mach number, the integrand varies as M e ®~ 5 . Since B~ 2.2 
for cases with a low wall temperature, the integrand decreases more rapidly with 
increasing at high Mach number, than in the case of a flat plate, and the ef- 
fects of an error in xg vanish more slowly. 


5 . 11 Coupling with the Inviscid Flow 

The method by which the boundary layer is coupled with the inviscid flow in 
NATA has been explained in Section 4. Briefly,' the boundary layer solution is 
generated step by step along with the inviscid solution. The inviscid flow vari- 
ables at each point are used in determining the boundary layer parameters n , d , 

5 , etc. , and the displacement thickness 8 at each point is used in the inviscid 
solution to provide a relation between the geometric area ratio A g and the effec- 
tive area ratio A e (Section 4.4). 

A fundamental problem in the calculation of coupled solutions for a confined 
inviscid flow and the boundary layer on the confining wall is that the system of 
equations for the coupled flow is unstable at supersonic Mach numbers (ref. 27, 
p. 377). If the boundary layer is computed using the Cohen-Reshotko integral 
method, as in the NATA code, this instability arises as follows: Assume that 

there is a small initial error in the Mach number gradient, dMe/df. If this 
gradient is too high, then the correlation parameter n (174) has too large a nega- 
tive value. Then the exponential factor e4-6n i n the correlation (189) for the 
incompressible form factor Hj nc is too small, and as a result, (1 + Hj nc ) is too 
small. Correspondingly, the displacement thickness 8 calculated from (188) is 
too low. When this erroneous 8 is used in the inviscid flow calculation, the 
effective area ratio computed is too high, so that the corresponding Mach number 
is too high. Then the next-computed value of the Mach number gradient, dM e /d£, 
is too high. Under some conditions, the error in the new dMg/df can be larger 
than the previous error which led to it. In such cases, the error is amplified 
from step to step, and the solution becomes unstable. Such instabilities were 
actually observed in some solutions run during the development of the NATA code. 

To avoid such instabilities, NATA uses a computational artifice. The incom- 
pressible form factor (189) is calculated from a smoothed value n of the correlation 
parameter n , rather than from the actual current value of n . In some earlier ver- 
sions of the code, n was defined as an unweighted average of n over the entire 
portion of the boundary layer upstream of the current flow point: 
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However, this has been found to be an unnecessarily severe approximation. In the 
current version of the code, ff is defined as an exponentially weighted average over 
the upstream part of the solution: 



Here a represents a characteristic averaging distance, such that n depends mainly 
upon the values of n at points lying between x-a and x . The coefficient b is 
chosen such that, for n = constant, (215) gives n = n . For x-xQ>>a , the value 
b = 1/a gives accurate results. The algorithm actually used in the code is de- 
rived from (215) as follows: Let tip be the n at the previous point of the solu- 

tion, and let Ax denote the current step Size. Then application of (215) to the 
current solution point gives, with b = 1/a f 



The integral on the right is approximated by treating n as constant over the step. 
Then 

n = ffp + n (1 — e — Ax/a) (217) 

This relation gives n as a weighted average of the n for the previous step and 
the n for the current step. 

The step size in NATA solutions is typically small in the throat region and 
quite large in the downstream region of high Mach numbers. An a value of a few 
times the local step size is required for stability. If a is many times the local 
step size, the solution remains stable but is affected to an inordinate extent by 
the use of n in place of n . Therefore, the averaging distance used in NATA is 
varied according to the following rule: 

a = R 0 Ag/w (218) 

where Rq denotes the throat radius, A g the geometric area ratio, and w a constant. 
Over a large part of the supersonic region, A^/Ax is approximately constant for 
frozen NATA inviscid solutions. In this region, (218) gives a/Ax ~ constant. 

Averaging of n over a fixed number of previous points in the solution would 
be simpler than (218) , but would cause the solution to depend upon the step size 
used in running the problem. The averaging prescribed by (217) and (218) makes 
the solution approximately independent of step size. 
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The smoothed value n of the correlation parameter is used only in calculating 
the form factor (and thus the displacement thickness). The shear parameter i and 
Reynolds analogy factor Ra are computed from the actual current value of n . Thus , 
the shear stress and heat flux are unaffected by the smoothing procedure, apart 
from any resulting errors in the inviscid solution itself. 


5.12 Geometric Options 

The basic geometric options available in NATA (two-dimensional nozzle, axi- 
symmetric nozzle, and rectangular channel) have been described and discussed in 
Section 4. The present section deals with the treatment of these options in the 
boundary layer calculations. 

The application of the Cohen-Reshotko integral method to two-dimensional and 
axisymmetric boundary layers has been discussed explicitly in preceding parts of 
the present section. The entire difference between these two cases is contained 
in the exponent j appearing in equation (171). If n' is calculated from equations 
(171), (172) and (173), and 6 from (185), then all of the boundary layer proper- 
ties are obtained simply by setting j = 0 for the two-dimensional case and j = 1 
for the. axisymmetric case. 

However, the third option (the rectangular channel) requires some further 
discussion. A channel has two pairs of nominally identical faces. If a channel 
face is of constant width, the boundary layer on it can be approximated as two- 
dimensional. If a face has variable width, the streamlines in the boundary layer 
converge or diverge. The effects of such a flow geometry upon the continuity 
equation for the -boundary layer can be taken into account by treating the boundary 
layer as axisymmetric in this case, but omitting the transverse curvature correc- 
tion. For example, a channel face which is widening in the downstream direction 
can be regarded as equivalent to the surface of a diverging axisymmetric nozzle 
whose circumference is equal to the width of the actual face at each axial position. 

Since the widths of the two sets of faces of a channel vary differently with 
axial position, two independent boundary layer calculations are performed in 
channel flow solutions. Both boundary layers are treated formally as if they 
were axisymmetric, i.e., using equation (170) to calculate n , For each layer, the 
equivalent radius r in (170) is taken to be the half-width of the channel face 
upon which the layer lies. No attempt is made to correct for the interaction of 
the boundary layers on adjoining channel faces in the comer regions. 


5.13 Examples and Discussion 

The present section exhibits selected results of NATA flow solutions includ- 
ing the boundary layer. These examples are intended to illustrate certain aspects 
of the NATA boundary layer calculations and to provide, by comparisons with exper- 
imental data, some indication as to the accuracy of the results. 

Figures 16 through 20 show some of the results of a non-equilibrium solution 
of the flow in a rectangular channel. The channel has a 2.54 x 5.08 cm (1x2 inch) 
cross section at the throat. In the downstream region, the 5.08 cm dimension is 
constant, and the other dimension increases from 2.54 to 53 cm at the exit, 145 cm 
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Figure 16 HEAT FLUX AND SHEAR STRESS FOR THE CHANNEL TEST CASE 
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Figure 17 PRESSURE AND MACH NUMBER FOR THE CHANNEL TEST CASE 
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Figure 1 8 GEOMETRIC AND EFFECTIVE AREA RATIOS FOR THE CHANNEL TEST CASE 
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Figure 19 MOMENTUM AND DISPLACEMENT THICKNESSES FOR THE CHANNEL TEST CASE 
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Figure 20 BOUNDARY LAYER CORRELATION PARAMETERS FOR THE CHANNEL TEST CASE 
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beyond the throat. In the region between 103 and 143 cm downstream of the throat, 
one of the wide faces of the channel is instrumented with pressure taps and flush- 
mounted heat transfer gages. 

The flow solution illustrated in figures 16 through 20 simulates a test with 
an actual channel at the NASA Johnson Space Center ARMSEF facility. The heated 
air flow was supplied by the Dual-Constrictor Arc. The conditions of the test 
were as follows: 

Mass flow = 0.0454 kg/sec (0.1 lb /sec) 

Plenum pressure = 0.919 atm 

Reservoir temperature = 6879° K 

Stagnation enthalpy = 25.1 MJ/kg (10799 Btu/lb) 

Reservoir molecular weight = 18.36 gm/mole 

The mass flow and plenum pressure are measurements. The other values were com- 
puted by NATA from these data. 

Figure 16 shows the variation of heat flux and shear stress with the axial 
coordinate,* , for this case, for the region from 0.1 m beyond the throat to the 
end of the channel. The curve labeled q w l represents the calculated heat flux to 
one of the broad faces of the channel. The three vertical bars represent heat 
flux measurements. Experimental data were available from three heat transfer 
gages at each of three axial positions. The ranges of values indicated by the 
vertical bars show the lateral variation of heat flux at each axial station. In 
most cases, the heat flux measured on the centerline of the channel face was higher 
than the values obtained from the off-axis gages. The main source of the varia- 
tion in heat flux was lateral non-uniformities in the flow. 

The curve labeled q W 2 represents the calculated heat flux to the narrow 
(5.08 cm wide) face of the channel. q w 2 is smaller than q w i because the boundary 
layer on the narrow face is thicker than that on the wide face. The curves 
labeled and r w 2 represent the calculated shear stresses on the wide and narrow 
faces of the channel, respectively. No experimental measurements of q w 2 » r w i , or 
1^2 are available. 

Figure 17 shows the variation of the static pressure ' the Mach number with 
axial position for the same case. The three circles rep - experimental 
measurements of the pressure on the centerline of one of i_ne wide channel faces. 

The Mach number at the channel exit is about 4.5. 

Figure 18 shows the geometric area ratio A_ and the effective area ratio A e 
as functions of axial position. At the exit, A e i s about 30 percent smaller than 
Ag. Figure 19 shows the momentum and displacement thicknesses. The momentum 
thickness 6\ on the broad face is smaller than that (# 2 ) on the narrow face of 
the channel because of flow-divergence in the boundary layer on the broad face. 

For x less than 0.155 meter, both displacement thicknesses, Sj and S 2 are nega- 
tive. In this region, the effective* cross sectional area of the flow (A^ ) is 
larger than the geometric cross section ( A^ ) . However, as shown in figure 18, 
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the effective area ratio A e is smaller than the geometric area ratio A g . These 
data are mutually consistent because the effective flow cross section A ' m at the 
sonic point is larger than the geometric cross section A'q at the throat. (See 
equations (113) and (114) . ) 6 

Figure 20 shows the correlation parameters-, nj and n 2 , for the boundary 
layers on the broad and narrow channel faces, respectively, together with the 
corresponding smoothed values and n 2 , calculated from equation (217). The 
smoothed values track the local values quite well. The local values nj and n 2 are 
shown for every step of the non-equilibrium solution in the region covered by the 
graph. Irregularities in step size between x = 1 m and the exit are due to in- 
sertion of steps at specified channel widths corresponding to the locations of 
heat transfer gages and pressure taps. 

Figure 21 presents a comparison of NATA heat flux calculations with experi- 
mental data on heat flux to the channel wall for eleven cases. The mass flows 
for these cases ranged from 0.028 to 0.091 kg/sec, the plenum pressures from 
0.544 to 1.985 atm, and the stagnation enthalpies from 10 to 47 MJ/kg. The verti- 
cal bars indicate the range of heat flux values observed at different lateral 
positions at each axial station. The dashed line is the locus of agreement be- 
tween the calculations and the experimental data. The unbroken line is a fit by 
eye to all of the data bars shown. Based on the unbroken line, the heat flux as 
given by NATA is too low by 28 percent, on the average. 

Figure 22 is a similar comparison of NATA static pressure calculations with 
experimental measurements using the centerline pressure taps in the channel. The 
flow cases are the same as in figure 22. Based on the curvefit provided by the 
unbroken line, the static pressures given by NATA are too low by 20 percent, on 
the average. 

One possible explanation of this discrepancy in static pressure would be 
that the effective area ratio, as computed in NATA, is too large because of errors 
in the boundary layer displacement thicknesses. To test this hypothesis, the 
ratio of experimental and calculated pressures is plotted versus the calculated 
displacement thickness on the broad face in figure 23 for a single station in the 
channel. If the errors in pressure were due to a systematic error in displacement 
thickness, this ratio would be close to 1 for small Sj , and would diverge from 
unity as 8^* increased. Figure 23, however, shows that the ratio has no significant 
correlation with displacement thickness. Thus, the hypothesized explanation is 
rejected. The apparent systematic error in static pressure is probably a result 
of lateral non-uniformities in the flow. The pressure data were all taken on the 
centerline of the channel face. If the off-axis pressures were lower than the 
centerline pressure, then the average pressure would be more nearly in agreement 
with the results of the NATA calculations, which are based on a uniform-flow model. 

The boundary layer calculations in NATA involve many approximations, includ- 
ing those listed below: 

(1) The assumption of uniform inviscid flow 

(2) The correlation assumptions in Cohen and Reshotko's method 

(3) The analytical curvefits (164), (165), (189), (194), (200), (201) to the 
results of similar solutions 
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Figure 21 COMPARISON OF EXPERIMENTAL AND CALCULATED VALUES OF HEAT FLUX 

TO THE CHANNEL WALL 
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Figure 22 COMPARISON OF EXPERIMENTAL AND CALCULATED VALUES OF 

STATIC PRESSURE 
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Figure 23 RATIO OF EXPERIMENTAL AND CALCULATED STATIC PRESSURE AS A FUNCTION 
OF CALCULATED DISPLACEMENT THICKNESS ON THE BROAD CHANNEL FACE 





(4) The assumed initial condition (Section 5.10) 

(5) The smoothing (217) of the correlation parameter to stabilize the coupled 

problem of the inviscid flow and the boundary layer 

(6) The linearity of the relation (164) assumed between the momentum param- 
eter N and the correlation parameter n 

(7) The neglect of the dependence of N upon the hypersonic parameter 

(8) The use of a Reynolds analogy factor, based on similar solutions, to 

calculate the heat flux 

(9) Approximations to gas properties used in the analysis, i.e., the perfect 
gas law and the viscosity-temperature relation (143). 


The combined effects of these approximations could account for sizeable errors in 
the results of 'the boundary layer calculations . In the context of these approxi- 
mations, the discrepancies actually seen in figures 21 and 22 appear moderate in 
magnitude . 
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6. EQUILIBRIUM AND FROZEN FLOW 


NATA performs two types of thermo chemical equilibrium calculation: 

(1) Equilibrium at specified temperature and pressure. This calculation is 
done by subroutine EQCALC, It is used in -the determination of the res- 
ervoir conditions and in the computation of conditions behind an equilib- 
rium normal shock on a model. 

(2) Equilibrium at specified temperature and entropy. This calculation is 
done by subroutine NEWRAP and is used in generating the solution for 
equilibrium flow of the gas through the nozzle. 

The technique used in the thermochemical equilibrium calculations at specified 
temperature and pressure is explained in Section 6.1. Section 6.2 discusses the 
inviscid equilibrium flow calculation. Section 6.3 discusses the coupling of 
the inviscid flow with the boundary layer on the nozzle or channel wall. Sec- 
tion 6.4 discusses the frozen flow solution, which is generated using techniques 
similar to those used in the equilibrium solution, with the constraint that the 
species mole fractions are constant at their reservoir values. Finally, Section 
6.5 explains the three options available in NATA for specifying the reservoir 
conditions, based on input of reservoir temperature and pressure, of reservoir 
pressure and total mass flow, and of stagnation enthalpy and total mass flow. 


6. 1 Thermoehemical Equilibrium at Specified Temperature and Pressure 
A chemical reaction can be specified symbolically in the form 

reactants products 


where the Mj, Mj represent chemical species and the Vj , v-' are the stoichiomet- 
ric coefficients. The condition for equilibrium of the reaction is (ref. 2, pp. 
284-285 and 952-953) 


reactants products 


( 220 ) 


in which pj , represent the chemical potentials of the reactant and product 
species. From equation (21) and the relation 


Pj = p Xj 


( 221 ) 
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between partial pressures pj , mole fractions Xi, and the gas pressure p, the 
chemical potentials can be expressed in the form 

f*i = (*i° + P +^ x i (222) 

in which the chemical potential ptj° at the standard pressure is a function only of 
temperature. 

Application of the equilibrium condition (220) to the reactions (5) for 
forming the dependent species (i - c + 1, .... n) from the independent species 
(m = l, C ) gives 


c 



m = 1 


(i = c + 1, . . . , n) 


(223) 


Substitution of equation (222) for the chemical potentials appearing in this 
relation gives an explicit expression for the mole fraction Xj of the i* depend- 
ent species in equilibrium with the independent species: 


c 



k = 1 


(224) 


where 


K] = exp 



(225) 


is the equilibrium constant. 

Conservation of the chemical elements can be formulated using equation (16): 


n 

x i- q i- X) x i ["<-=. i - "i 'fi-.- 1 '] 

i =c + l 


(226) 


Here the qj are the normalized coefficients (10) expressing the elemental composi- 
tion of the gas in terms of the independent species, and the Xj (for j = 1, — , c) 
are the mole fractions of those species. The (n-c) equations (224) together with 
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the c equations (226) provide a system of n equations for determining the mole 
fractions X m ( m = 1,..., n) 0 f all the species. However, if equations (224) are 
used to eliminate the mole fractions Xj of the dependent species from (226) , 
equations (226) become a system of c equations for the mole fractions Xj of the 
c independent species. This formulation of the thermochemical equilibrium prob- 
lem thus reduces the number of equations to be solved from n (the number of 
species) to c (the number of elements). Since the right-hand side of (224) con- 
tains a product of mole fractions X m raised to various powers vj _ c m , the system 
of equations is, in general, non-linear. 

In NATA, these equations for the mole fractions of the independent species 
are solved using the Newton- Raphson method. Equations (226) are rewritten in 
the form 


F: (X V 


XJ = 


qj- 


U 

Xi ['i-M-’i <1 '?-c- 1 >] - 


X; 


i = c + 1 


(227) 


The component mole fractions in the r th iteration will be denoted by Xj r . Then, 
to first order, 


VV* 1 ' - f, o k ') * 



<V +I -V> 


(228) 


where (dFj/dX,,,) 1 denotes the derivative of Fj with respect to X m , evaluated using 
the mole-fraction values of the rd> iteration. This derivative can be determined 
analytically by differentiating (227) : 



i = c + l 


(229) 


since the qj , vj_ c : , and v*. _ c are all independent of the mole fractions for a 
given gas model. The subscripts j and m refer to independent species, while i 
refers to a dependent species. Thus 


dXj _ ( 1 

dx m S i m ~ ) 
m / o 


for 

j =m 

for 

j m 


( 230 ) 
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and from logarithmic differentiation of (224f, 


«?X; 


ax, 


k = l 


v i - c, k %m v i — c,i 


X, 


( 231 ) 


Substitution of these expressions into (229) gives 


ax„ 


u 

. x- 


8 - 


jm 


(232) 


i=c + l 


The object of the Newton-Raphson iteration is to make the Fj (X k ) all equa.l 
to zero. Hence, in equations (228), the Fj (X k r + 1 ) are set equal to zero, for 
j = 1 to c. This gives the following system of equations: 


c 




Fj <X k r ) 


(233a) 


where the h„ r are the relative corrections to the mole fractions: 


h m r s 
m 


v r+1 Y 1 
A m A m 


X„ 


such that the corrected mole fractions are given by 

X r + 1 = X „ r (1 +h r ) 
m m ' m ' 

(233a) is a system of c linear, inhomogeneous equations for the h m 
of (232) into (233a) gives the matrix of coefficients in the form 


(234) 

(235) 

Substitution 


- X„ 



u 


v i - c, m + ^jm 


i = c + l 


(233b) 


In the special case of a gas model in which there are no dependent species 
(i.e., in which the number of species n is equal to the number of chemical ele- 
ments c ) , the above formulas remain valid if the sums running over the dependent 
species are eliminated. These are the sums in which the summation index runs 
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from c + 1 ton. Equation (226) shows that, in this special case, the mole 
fractions Xj are equal to 9j and thus are constant. Equations (227), (233a), 
and (233b), with the sums over dependent species eliminated, show that the 
Newton-Raphson procedure gives this solution for the independent mole fractions 
correctly on the first iteration. 

In NATA, the system of equations (233) is solved for the corrections h m r , 
at each stage of the iteration, by calling a subroutine (DSMS0L) for simultaneous 
solution of a system of linear equations. The iteration is continued until the 
h m r are all less than or equal to 10-6 in absolute magnitude. Once the iteration 
has converged, the entropy, enthalpy, and mean molecular weight of the equilib- 
rium gas mixture are computed from equations (30b) , (26) , and 



(236) 


The gas density is then computed from the ideal gas relation 


p W 

V? 


(237) 


Finally, the density and enthalpy are corrected for gas imperfections using 
equations (79b) and (81) . 


It should be noted that the mixture enthalpy given by (26) , 


H = 



X i H i 


(238) 


is the enthalpy per mole of mixture. The specific enthalpy (enthalpy per unit 
mass) is 



(239) 


where W is the mean molecular weight (236). In the internal calculations of 
NATA, the species molar enthalpies Hj are computed in the non-dimensional form 

H j 

SHJ (J) = (240) 

R 0 T 0 

where Rq is the universal gas constant and T 0 the reservoir temperature, and 
the specific enthalpy is given as 
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CH 


(241) 


W 0 h 

R o T o 


where Wg is the molecular weight in the reservoir. From equations (238) to 
(241), 


w o v"* 

CH = — y Xj . SHJ (J) (242) 

J = 1 


6.2 Inviscid Equilibrium Flow 


The governing equations for steady, one-dimensional adiabatic inviscid flow 

are 


P u A e = P* u * A e* = const ' 


(243) 


u du + 


dp 


0 


h 



const. 


(244) 

(245) 


where p denotes the density, u the flow velocity. A' the cross sectional area, 
p the pressure, and h the specific enthalpy. In the continuity equation (243), 
the asterisks indicate sonic-point values. The Euler equation (244) assumes that 
viscous stresses and body forces (such as electromagnetic and gravitational 
forces) are negligible. The energy equation (245) assumes that the flow is 
adiabatic. Subject to these assumptions, these equations are applicable to any 
steady, one-dimensional flow whether the gas is chemically frozen, in equilib- 
rium, or undergoing reactions with finite rates. 

Equations (244) and (245) can be combined to give 



9 

From the first law of thermodynamics, 

Tds > de + p d 



(246) 


(247) 
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and the definition of enthalpy 


h 


P 

e + — 


P 


it may be shown that 


dp 

Tds > dh — 

P 


( 248 ) 


(249) 


Equation (246) thus implies that 

Tds > 0 (250) 

i.e., the specific entropy of the gas either remains constant or increases in a 
flow obeying equations (244) and (245). The condition for ds = 0 is that the 
flow be reversible in the thermodynamic sense. This means that if the flow 
retraces its sequence of velocities (without passing through a shock), the gas 
is restored to its initial state. The flow id reversible, in this sense, if 
the gas obeys a rate- independent equation of state, p * p (p, h), such that the 
density is uniquely determined when the pressure and enthalpy are specified, A 
gas which is everywhere in local thermo chemical equilibrium obeys a constitutive 
relation of this type. A gas whose chemical composition is frozen also obeys a 
rate- independent equation of state. Thus, both equilibrium and frozen flows 
governed by equations (243) to (245) are isentropic (ds= 0 ). On the other hand, 
non-equilibrium flows with finite reaction rates are Irreversible, because the 
gas composition depends not only upon the pressure and enthalpy but also upon 
the past history of the gas sample. Thus, from (250), the entropy increases in 
such non-equilibrium flows. 

All of the flow solutions calculated by NATA are assumed to start from a 
state of thermochemical equilibrium in an upstream reservoir. The composition 
and state of the gas in the reservoir are computed using the method described 
in Section 6.1. If the reservoir temperature and pressure are input-specified, 
this method is applied directly. If other options for input specification of 
the reservoir conditions are employed, the reservoir temperature and pressure 
are determined using the iterative techniques explained in Section 6.5, and the 
method of Section 6.1 is again used to compute the gas composition and state in 
the reservoir. 

The flow conditions that can be reached by an equilibrium expansion from a 
specified reservoir condition form a one-parameter family of states. In NATA, 
the gas temperature is taken to be the independent variable. A sequence of 
equilibrium flow states is generated by decrementing the temperature, starting 
from the reservoir value. At each temperature, the species mole fractions and 
the static pressure are determined from a thermochemical equilibrium calculation 
with the supplementary condition that the gas specific entropy s be equal to its 
value s 0 in the reservoir. The conditions for thermochemical equilibrium are 
equations (224) and (227). The specific entropy may be obtained by dividing 
the molar entropy (30b) by the local molecular weight w: 
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( 251 ) 
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Thus, the condition s = s 0 can be written, with the aid of (236), in the form 


I 


S-° s, 


j = 1 


0 ^0 


'0 w j 

- - fcx X: 

Rn » 


£ri p — 0 


(252) 


in which Xj iA the mole fraction of the j th species, Wj the species molecular 
weight, Sj° the species molar entropy at standard pressure, and R 0 the universal 
gas constant. Equations (227) and (252) constitute a system of c + 1 equations 
to be solved for the equilibrium mole fractions and the static pressure p at the 
specified temperature and specific entropy. As in Section 6.1, these equations 
are solved by the Newton-Raphson method. The solution is carried out by subrou- 
tine NEWRAP. The additional derivatives required by incorporation of the addi- 
tional equation (252) and the additional unknown P are as follows, for j = 1, .... c 
and m = 1 c : 



dp 


53 _q ’ (v *i-c- 1) ]“7' (v *i-c- 1} 


i= c+ 1 


(253) 


aF c + l V 

= R 0 


Rn 


- ^ X m - 1 



i = c + 1 


(254) 



s 0 W i 



( 255 ) 
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Equation (253) is derived by differentiating equation (227) with respect to the 
pressure, p . For the independent species (j) , 



(256) 


since the mole fractions Xj (j = 1, .... c) and the pressure p are all independent 
variables. For the dependent species (i) , from (224), 


dx i 

dp 



(257) 


Equation (254) is obtained by differentiating (252) with respect to the mole 
fraction X m of an independent species, and using equations (230) and (231). 
Equation (255) results from differentiation of (252) with respect to p , when 
(256) and (257) are taken into account. 

The system of Newton-Raphson equations for the state of the gas at a point 
in an equilibrium flow may be written 
c + 1 

-£ v 

m = 1 

where j and k run from i to c + 1 , 

Y m r - XJ (“ = !. •••• c > (259a) 


Ifi 

Jy„ 


h m f = F j (Yk> 


(258) 


W = P r 


and 

Y t + 1' y * 



(259b) 


(260) 


For m = 1 to c and j = 1 to c , the matrix of coefficients in (258) is given by 
(233b) . For m = c + 1 and j =1 to c, 



0- C)i - ij x >] x i (y\_ c - l) 


(261a) 


- 107 - 



For m = 1 to c and j = c + 1 , from (254) 
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(261b) 


+ | tn Xj ) :(v j _ c - 1) (261c) 


i = c + l 


At each stage of the Newton-Raphson iteration, the system of c + 1 linear 
equations (258) is solved for the relative corrections hnf by calling subroutine 
DSMS0L. The improved estimates for the independent species mole fractions X n and 
the pressure are then obtained from 


v r + 1 
m 


Y r 

m 


« +h m r ) 


(m = 1, c + 1) 


(262) 


The iteration is continued until all of the are smaller than or equal to 

1(T 6 . 


Once the equilibrium mole fractions and the static pressure have thus been 
determined, the molecular weight, density, and specific enthalpy are computed 
from equations (236) , (237) , and (242) . The density and enthalpy are corrected 
for effects of gas imperfections using equations (79b) and (81) . The flow 
velocity is then determined from equation (245) : 

u = / 2(h 0 _ h) C 263 ) 

and the mass flux 

m = pu (264) 


is computed. 

Up to this point, the calculations have proceeded without reference to the 
nozzle geometry or to the location in the nozzle at which the computed flow con- 
ditions occur. The flow conditions are related to the nozzle geometry only 
through the continuity equation (243) . This equation can be solved for the 

area ratio A : 
e 

a ' P u m 

A = = J! 1 = (265) 

e A' pu m 

e * 
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Thus, if the sonic mass flux m, = P* u t is known, the area ratio can be calculated 
from the local mass flux (264) . The location of the flow point within the nozzle 
can then be determined by inverting the area ratio relation A e = f(x) , based on 
the nozzle or channel geometry (Section 4) , to obtain the axial position x cor- 
responding to the calculated area ratio. 

Since the area ratio A e has a minimum at the sonic point, the sonic mass 
flux m, is the maximum mass flux occurring anywhere in the flow. It is determined 
from a preliminary calculation (performed in subroutine NRMAX) , before the main 
equilibrium flow solution. In NEMAX, the temperature is repeatedly decremented, 
and the equilibrium mass flux is computed at each stage, until m passes its 
maximum value. An iteration is then carried out to determine the temperature 
and other flow conditions corresponding to the maximum mass flux. The resulting 
m* value is then used in the determination of the area ratio from (265) during 
the main equilibrium flow solution. 


6.3 Boundary Layer Effects 

Inviscid flow is, of course, only an idealization of the behavior of real 
gases. In an actual wind tunnel, there is a boundary layer on the nozzle wall, 
within which viscous stresses and heat transfer are important. Upstream of the 
throat and in the throat region, this boundary layer is thin and has small (but 
not negligible) effects upon the free stream flow. Downstream, in the region 
of high Mach number and low density, the boundary layer thickness can become 
comparable with the nozzle radius or the gap between channel faces. In such 
cases, the effects of the boundary layer upon the free-stream flow can become 
quite important. 

Depending on the flow conditions, the boundary layer can be either laminar 
or turbulent. Within the operating envelopes of existing NASA/ JSC arc-heated 
wind tunnels, the layer is usually laminar. NATA contains an approximate calcu- 
lation of laminar boundary layer development on the nozzle wall, based on the 
Cohen-Reshotko integral method (Section 5) . The boundary layer properties , 
including the displacement thickness 8 * , are computed step-by-step, beginning 
at an upstream starting point, along with the inviscid flow solution. The 
displacement thickness alters the effective geometry of the nozzle or channel.* 
If it is positive, it reduces the effective cross sectional area of the flow. 

If it is negative, as it normally is upstream and in the throat region, it 
increases the effective area. 

Relations between the effective area ratio (allowing for the displacement 
thickness) and the geometric area ratio have been given in Section 4. These 
relations involve the displacement thickness '8 at the nozzle throat. This 
quantity is not known until the inviscid flow solution and the boundary layer 
calculation have reached the throat. Thus, the solution in the region upstream 
of the throat must be generated without accurate knowledge of the relation 
between the effective and geometric area ratios in this region. Several dif- 
ferent techniques for circumventing this difficulty have been tried during the 


# ln the case of a channel, the flow is confined by two pairs of walls with different geometry. In this case, NATA performs two 
separate boundary layer calculations, and there are two displacement thicknesses at each flow point. 


- 109 - 



development of NAT A. Most of these techniques involved the use of an approximate 
value for the displacement thickness at the throat, based upon some type of pre- 
liminary solution up to the throat. Techniques of this type involved relatively 
complex programming and, in some cases, led to unreliable code operation in the 
throat region, where the solution is extremely sensitive to small errors in the 
area ratio. 

The technique used in the current version of NATA is simple and reliable, 
and gives results of acceptable accuracy. Upstream of the throat, the relation 
between the effective area ratio and the axial coordinate x is approximated by 
neglecting the difference between the effective and geometric area ratios in 
this region. Since the boundary layer is thin in comparison with the nozzle 
diameter upstream and near the throat , this approximation is equivalent to a 
small change in the nozzle shape in the upstream region. Such a change has 
little effect on the displacement thickness at the throat. Downstream of the 
throat, in the region where the boundary layer thickness can become large, the 
relation betwden the geometric and effective area ratios is calculated using 
the local displacement thickness, so that the boundary layer and the inviscid 
flow are coupled. 

The effect of the displacement thickness upon calculations of reservoir 
conditions from the total mass flow using sonic flow analysis is discussed 
below, in Section 6.5. This effect is not negligible, and is taken into account 
by an iteration. 

The preceding discussion of boundary layer effects is phrased in general 
terms, and is applicable to all three types of flow solution. In the specific 
case of equilibrium free-stream flow, the solution including boundary layer 
effects is generated in the following way: 

(1) As in the purely inviscid case discussed in Section 6.2, the gas 
temperature is taken as the independent flow variable. The tempera- 
ture is repeatedly decremented, starting from the reservoir, and at 
each stage the composition, thermodynamic properties, and flow 
velocity are determined by a Newton-Raphson solution of equations 
(227) and (252). 

(2) Upstream of the throat (i.e., for temperatures above the sonic tempera- 
ture T„ ) , the axial coordinate corresponding to each flow point is 
determined from equation (265) by assuming that the effective area 
ratio A e is equal to the geometric area ratio Ag, That is, x is 
determined by solving the equation 


A g (x) 



(266) 


with the aid of subroutine FINDX. 

(3) For each flow point, after x has thus been determined, the boundary 
layer calculation for the point is carried out as described in 
Section 5. 


- 110 - 



(4) This procedure is continued until the throat (x = 0) is reached. The 
displacement thickness at the throat (8** ) is thus determined 
approximately . 

(5) Beyond the throat, the flow conditions are again determined by decrement- 
ing the temperature and solving equations (227) and (252) , but the geo- 
metric area ratio A g is obtained from equation (127), (131), or (135), 
depending on the type of nozzle geometry. The axial coordinate x is 
then determined by inverting the geometric area ratio relation A g (x) for 
the nozzle or channel. These operations are performed in subroutine 
AGS0LN. Since the formulas relating A g and A e involve the displacement 
thickness 8*, and since 8* is obtained from the boundary layer calcula- 
tion which requires the coordinate x, the computation of Ag is based on 
the effective area ratio A e for the current flow point and the displace- 
ment thickness 8* at the preceding flow point. After the boundary layer 
calculation for the flow point has been performed, subroutine AGS0LN is 
called again with the effective area ratio and displacement thickness 
for the current flow point. The resulting Improved value of x is the 
one printed by the code. The error resulting from use of this approxi- 
mation is small so long as the change in 8* in one step of the flow 
calculation is sufficiently small in comparison with the nozzle 
diameter . 


6.4 Frozen Flow 

Equilibrium flow is the limiting case of the flow of a reacting gas mixture 
in which the reaction rates are infinitely high. Frozen flow is the opposite 
limit, in which the reaction rates are zero. Thus, in frozen flow, the species 
mole fractions are assumed to be constant and equal to their values in the up- 
stream reservoir. 


As in the case of equilibrium flow, the frozen flow solution is generated 
in NATA by decrementing the temperature, starting from the reservoir value. At 
each temperature, the flow conditions are computed as follows: The enthalpy is 

computed from (242) . The static pressure is determined from the condition that 
the flow be isentropic; from (251), this condition can be formulated 


P 





(267) 


in which s 0 denotes the specific entropy in the reservoir and W 0 the molecular 
weight, which is everywhere equal to the reservoir value since the flow is 
frozen. Once the pressure is known, the perfect-gas density can be calculated 
from (237) . For a frozen flow, (237) becomes 


P_ = P_ 

P0 Po T 


(268) 
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The density is calculated from (268), and the density and enthalpy are then cor- 
rected for the effects of gas imperfections using equations (79b) and (81). 

These calculations are all performed, in subroutine PR0P. 

As in the equilibrium case, the inviscid flow conditions are all computed 
without reference to the nozzle geometry. The sonic mass flux for frozen flow 
is determined by a preliminary calculation, performed before the main frozen 
flow solution, in which the temperature is adjusted repeatedly to find the value 
corresponding to the maximum mass flux. During the main frozen solution, at 
each temperature, the flow conditions are computed as described above and the 
effective area ratio A e is then calculated from the continuity equation in the 
form (265) . The location in the nozzle at which the calculated flow conditions 
occur is then determined from the effective area ratio. If a solution neglecting 
the boundary layer is being generated, the geometric area ratio Ag is equal to 
A e , and the axial coordinate x is obtained by inverting the area ratio relation 
Ag = Ag(x) for the nozzle. If the solution includes boundary layer effects, the 
technique for determining x from A e and the displacement thickness 8 * is exactly 
the same as in the case of the equilibrium solution, described in Section 6.3. 


6.5 Determination of Reservoir Conditions 


NATA provides three options for input specification of the reservoir condi- 
tions. In all three cases, it is assumed that the gas in the reservoir is in 
thermochemical equilibrium. In the. first option, the code user specifies the 
reservoir temperature and pressure directly. The composition and thermodynamic 
properties of the gas in the reservoir are then computed as described in 
Section 6.1. 

In the second option, the input quantities are the reservoir pressure and 
the total mass flow. The rationale supporting this choice of parameters is that 
the total mass flow is much easier to measure than the reservoir temperature. 

A meaningful reservoir pressure measurement can also be obtained easily in wind 
tunnels equipped with a plenum chamber within which a condition of flow stagna- 
tion is approximated. In this second option, the reservoir temperature is 
determined, by an iteration, based on the condition that the product of the 
effective throat area and the equilibrium sonic mass flux be equal to the input 
total mass flow. For each successive trial value of the reservoir temperature, 
the equilibrium reservoir conditions are calculated as explained in Section 6.1, 
and the equilibrium sonic mass flux m* is then computed as described in Section 
6.2. If boundary layer effects are to be neglected in the flow solution, the 
effective throat area is equal to the known geometric cross section of the 
nozzle at the throat. If boundary layer effects are to be included, the effec- 
tive throat area involves the displacement thickness at the throat, as shown 
in equations (125), (129), and (133). In this case, the reservoir temperature 
is determined initially neglecting boundary layer effects. A preliminary 
equilibrium flow solution from the reservoir to the throat is then computed to 
determine an approximation to the displacement thickness at the throat. The 
calculation of the reservoir temperature is then repeated, using an effective 
throat area which includes the effect of the boundary layer displacement 
thickness. The main flow solutions for frozen, equilibrium, and non-equilibrium 
flow are then carried out assuming this revised reservoir temperature. 
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The third option for specifying the reservoir conditions is based on input 
of the total mass flow and the stagnation enthalpy. The enthalpy can be deter- 
mined from energy-balance data, and is a more meaningful piece of information 
than the upstream pressure in wind tunnel configurations lacking a stagnation 
region between the arc heater and the throat, especially if the upstream flow 
is highly rotational. When this option is used, NATA determines the stagnation 
temperature and pressure in a fictitious upstream reservoir by a double itera- 
tion, based on the conditions that the reservoir enthalpy be equal to the input 
stagnation enthalpy, and that the product of the effective throat area and the 
equilibrium sonic mass flux be equal to the input total mass flow. If boundary 
layer effects are to be included in the flow solution, a preliminary calculation 
to determine the displacement thickness at the throat is performed prior to the 
main flow solutions, as in the second option. 

In the second and third options, the determination of the reservoir condi- 
tions is based in part on calculations of equilibrium flow from the reservoir 
to the throat, regardless of whether the main flow solutions to be generated 
are to be equilibrium, frozen, non-equilibrium, or some combination of these. 

The rationale for this procedure may be outlined as follows: 

(1) When two or three different types of flow solution are computed for the 
same problem, the user's purpose is to assess the importance of the re- 
actions by comparing solutions based on infinite, finite, and zero 
reaction rates. Such an assessment would be obscured if the different 
types of solution started from different reservoir conditions, as they 
would if separate reservoir calculations were performed for frozen, 
non-equilibrium, and equilibrium flow. 

(2) The equilibrium-flow reservoir calculations for the second and third 
options already consume a significant amount of computer time, ranging 
from a few seconds to over a minute per case, depending on the 

gas model and the option. Similar calculations assuming non-equilibrium 
flow would use a great deal more time, especially since the non- 
equilibrium solution is often forced, by stability requirements, to 
take very small steps in the upstream region. 

(3) In most of the cases to which NATA is applied, the non-equilibrium 
solution (which simulates the actual flow most closely) is approximated 
reasonably well by the equilibrium solution in the region upstream of 
the throat, because the flow starts from an equilibrium state in the 
reservoir and the pressure remains fairly high until the throat has 
been passed. 

The reservoir condition calculations for the second and third options are 
controlled by subroutine RESTMP, 
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7. NON-EQUILIBRIUM FLOW 


The chemical reaction rates in a real, high-temperature gas mixture are 
finite. Thus, equilibrium and frozen flow are only limiting cases, which under 
some circumstances may approximate the flow in portions of the nozzle. The 
actual flow in an arc-heated wind tunnel is a non-equilibrium expansion. The 
present section explains how non- equilibrium flow is calculated in the NATA 
code.* Section 7.1 presents the flow relations and reaction rate equations which 
are assumed to govern the flow, both for the normal case of chemical non-equilibrium 
and for the case in which the free electrons and bound electronic excited states 
are also assumed to be out of equilibrium with the heavy species. Section 7.2 out- 
lines the overall method of solution used in NATA. Special techniques used to 
calculate portions of the non-equilibrium flow are then described and analyzed in 
Sections 7.3 to 7.5. These techniques are an inverse method used in the upstream 
region, a perturbation method used to start the solution in the region just down- 
stream of the reservoir, where the flow is still nearly in equilibrium, and the 
numerical integration method used to compute the flow farther downstream, where 
some of the reactions are appreciably out of equilibrium. Finally, Section 7.6 
explains how the non-equilibrium inviscid flow is coupled with the boundary layer 
on the nozzle or channel wall. 


7. 1 Governing Equations 

The NATA code contains two different treatments of non-equilibrium flow. 

The first is a conventional single-temperature model in which, at each station 
in the nozzle, the kinetic translational temperatures and the excitation temper- 
atures of all species are assumed equal. In this model, only the species concen- 
trations are allowed to depart from equilibrium. The compiled-in air models and 
planetary atmosphere models are of this type. For brevity, these may be referred 
to as chemical non-equilibrium models. 

The second type of non-equilibrium treatment implemented in NATA is a two- 
temperature model. In this case, in addition to non-equilibrium of the species 
concentrations, the electronic degrees of freedom are assumed to be out of 
equilibrium with the translational and vibrational motions of the atoms, ions, 
and molecules. The velocity distributions of these heavy particles are assumed 
to be Maxwellian at a temperature T, while the translational temperature of the 
free electrons is allowed to have a different value, T e . The electronic excited 
states of some chemical species are treated as separate physical species, so that 
the populations of these states need not be in equilibrium at either of the 
temperatures, T or T e . The governing equations include terms representing energy 
transfer between the electrons and heavy particles. In addition, radiative losses 
from the plasma are taken into account. The compiled-in models for helium and 
argon ar._ cf this type. Such models will be referred to as electronic non- 
equilibrium models. 

The first part of this section presents and discusses the governing equations 
for the chemical non-equilibrium models. The second part deals with the elec- 
tronic non-equilibrium models. Part 3 explains the technique used to maintain 
the elemental composition of the gas mixture. 

•Available techniques for analyzing non-equilibrium nozzle flows have been reviewed by Bray (ref, 31) and by Hall and Treanor 
(ref. 32). 
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7.1.1 Chemical Non-equilibrium 


The gas flowing through the nozzle is assumed to consist of n chemical 
species which can undergo r chemical reactions. The reactions can be repre- 
sented in the form 



1j M ) 



( 269 ) 


In this formula, the Mj symbols represent the chemical species, and the matrices 
vjj , are the stoichiometric coefficients. The forward and reverse reaction 
rate constants are denoted by kfj , k r j , respectively. The gas is assumed to obey 
the equations of steady, one-dimensional, adiabatic inviscid flow,’ namely 
equations (243) through (245) . These equations may be written in differential 
form as 
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d^nu 
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(270) 
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(271) 
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du 


(272) 


Equations (270) to (272) are the continuity, momentum, and energy equations, 
respectively. The gas mixture is assumed to be ideal, so that the equation of 
state is 


P 


pR 0 T 

w - 


(273) 


in which the molecular weight W is given by equation (236) : 



(274) 


in terms of the mole fractions Xj and molecular weights Wj of the individual 
species. Also, the specific enthalpy h is given by 


h = 




(275) 


where Hj is the molar enthalpy of the j* species, a function of temperature. 
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Inspection of equations (270) through (275) reveals the presence of n +6 
dependent variables: p, u , P, h, T, W, and Xj for j = 1 to n, (The flow area 

A' e is assumed, for now, to be a known function, A' e (i) of the position coordinate 
x.) Since there are just six equations, it may be seen that n additional re- 
lations are required, to yield a determinate system. The required additional 
relations are the rate equations specifying the changes in the species mole 
fractions, Xj , or concentrations, yj . According to equation (61), the net rate 
of change of the molar concentration [M- ] of species j due to the i c ^ reaction 
is ’ 

n 

v ik TT 

- k ri [[ [M k l 
k = 1 


(276) 


I d [ Mj ] | 
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In this equation, [M k ] has units of moles/cm^. The forward and reverse rate 
constants, k fi and k ri , are assumed to be connected by the detailed balancing 
relation (62) , 



(277) 


in which the equilibrium constant Kj is given by equation (68) ! 


K i = 
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Here 
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(278) 


(279) 


(280) 


Some manipulation is required to bring equation (276) into a form suitable 
for use in the non-equilibrium flow calculations. First, in the NATA non- 
equilibrium solution, the amounts of the species are expressed in terms of the 
specific molar concentrations yj , which are related to the mole fractions by 
equation (1) and which have units of moles per gram of mixture. The volume con- 
centrations [Mj] are related to the y j by 

[Mj] = p yj ( 28 D 
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in which p is the gas density. Second, since NATA deals with steady flows, the 
time derivative in (276) must be converted into a space derivative using 


d d 

dc dx 


(282) 


Substitution of equations (279), (281) and (282) into (276) gives 


pa 


dx 
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(283) 


in which 



(284a) 


(284b) 


Equation (276) gives only the changes in [ Mj ] due to the chemical reaction, not 
those resulting from changes in the volume of the gas sample. Thus, the ld/dt}. 
operator on the left in (276) refers to changes in species concentrations at con- 
stant gas density, and it is correct to set 



as was done in the derivation of equation (283)*. 

The total rate of change of the specific concentration y ■ for the j* species 
may now be obtained by summing equation (283) over all of the reactions: 


pa 


dyj 

dx 




( 286 ) 


•See Section 7.1 .2 (below) for a more detailed derivation of (283). 
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With the aid of the detailed balancing relation (277), equation (286) can be re- 
written in the form used in KATA: 


r 
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(287) 


where 
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(288) 
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(289) 


The five flow equations (270) to (273) and (275) can be reduced to two by 
eliminating the derivatives du/dx , dp/dx, and dh/dx . Combination of the con- 
tinuity equation (270) with the energy equation (272) gives 


d Enp d?nA' e 1 dh 
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From (275), 
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Also, from (31), 
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Hence, (290) can be rewritten in the form 
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Differentiation of the equation of state (273) gives 
dp R 0 T dp R 0 p dT pR 0 T dW 

dx W dx W dx dx 

Substitution of (272) and (294) into the momentum equation (271) gives 

dh R 0 T dinp Rq dT RqT dW 
+ + — = 0 (295) 

dx W dx W dx d x 

The quantity dw/dx can be expressed in terms of the concentration gradients dyj/dx 
by differentiating the relation 



(296) 


which can be derived simply by summing equation (1) over the species. Differen- 
tiation of (296) gives 



(297) 


Substitution of (291), (292), (296), and (297) into (295) gives 

E d y: dT R 0 T d in P 

(R o T - H i> IT y i (r °- c pP IT + — — " 0 

j=i i-i 


(298) 


Addition of equations (293) and (298) gives, finally. 



(299) 


The basic system of differential equations for the chemical non-equilibrium model 
is (287), (293), and (299). This is a system of n +2equations for the dependent 
variables p, T, and yj for j = 1 to n. 


- 120 - 



In NATA, these equations' are expressed in finite difference form and solved 
numerically. The integration technique is discussed in Section 7.5. At each 
point x in the flow, the concentration derivatives dyj/dx are computed from 
equations (287), and the derivatives dT/dx and d p/dx are obtained by simul- 
taneous solution of equations (293) and (299). The conditions at a pointx+Ax 
are then calculated from the flow variables and their first-order derivatives at 
x, using a modified Runge-Kutta technique. Once T, p , and the yj have thus been 
determined at the new point, the specific enthalpy is computed from (275) and the 
velocity from the integral form of the energy equation (245) , To ensure accurate 
conservation of mass, the density is recalculated from the flow velocity and the 
area ratio at the current flow position x, using the continuity equation (265). 


The molecular weight is then computed from (296) , and the pressure from the 
equation of state (273). Also, the specific entropy is calculated from equation 


(251), and a Mach number is calculated based on a "speed of sound" given by 
[(dp/dx)/(dp/dx)l 1/2 . From (270) and (271), 



(300) 


The above defined speed of sound is not, in general, the speed of propagation of 
small disturbances in a relaxing medium (ref. 33). 


7.1.2 Electronic Non-equilibrium 

Measurements of the electron temperature T e in ionized nozzle flows generally 
show that the equilibrium between T e and the gas temperature breaks down at some 
point in the expansion. (See ref. 34, for example.) The importance of including 
this phenomenon in flow calculations depends upon the intended application of the 
results. Aerodynamic forces and heating are not likely to be much affected by 
electronic non-equilibrium, but the electron density, the populations of excited 
states, and the spectrally resolved radiation emitted by the gas all depend 
sensitively on the electron temperature. 

Equations for steady-state, quasi-one-dlmensional flow of a plasma with un- 
equal electron and heavy-particle temperatures have been formulated and solved 
in several previous studies (refs. 35 through 37). Since there appears to be 
some disagreement in the literature as to the correct form of such equations, 
the applicable relations are derived here by formulating conditions for mass, 
momentum, and energy conservation for each species. 

Mass Conservation . - Let rj represent the mass of the j th species produced 
per unit volume per unit time by reactions occurring in the flow. Then the mass 
conservation equation for species j may be written by considering the mass balance 
for a small volume element in the flow of thickness Ax. The total mass of the jth 
gas entering this volume element per unit time through the upstream boundary at 
x = x is then given by pj uj A' e evaluated at x, while the mass leaving the volume 
through the downstream boundary at x = x + Axis pj uj A' evaluated at x = x + Ax, 
where pj is the mass density of the j* species in tne gas, uj is its mean flow 
velocity, and is the cross-sectional area of the flow at the given value of x. 
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Since the net mass of species j leaving the volume element between x and x+Ax 
due to flow through the ends must be just equal to the mass of j produced in the 
volume element by the source term " t . , the mass balance for the volume element 
becomes ! 




x + Ax 


ij A g dx 


(301) 


Taking the limit of equation (301) as Ax-*-0 yields the mass conservation 
equation for the j* species in steady quasi-one-dimensional flow. 


d . 

— (pj Uj A e ) = rj A e 


(302) 


Momentum Conservation . - To derive the momentum conservation equation for 
the j® species, consider a test mass of the j th gas which is instantaneously 
located between the planes x and x + Ax and which is moving with the local flow 
velocity uj . The total momentum of this gas element is then 


P 


x + Ax 



X 


P j u j a; dx 


(303) 


while the force acting on the element is 

x+ Ax 


F = (p 


x + A, JA . 

- (p i 7 


fj a ; dx 


(304) 


Here, the first two terms on the right hand side of equation (304) are the net 
pressure force of the j* gas acting across the end planes of the volume element 
at x and x + Ax , the third term is the x-component of the force between the nozzle 
wall and the j* gas and the final term is the total force on the test mass of the 
j* gas due to interactions with other gaseous species, where fj represents the 
total body force per unit volume on the j th gas due to the interactions. Apply- 
ing the momentum conservation condition F = dP/dt to the test mass then yields 
the momentum balance for the element as 


- 122 - 



(305) 


(Pj - (p, 


/ 


x+ Ax 




p i ~ ^ + 


x+ Ax 

/ 


fj a; dx 


(p i A;U i ) X + A* - ( Pj A e u j 2) _ 


where the right-hand-side of equation C305) follows directly from equation (303) , 
since the limits of integration in (304) are assumed to move with the local flow 
velocity uj . Taking the limit of equation (305) as A x— ► 0 now yields the momentum 
conservation equation for the j th species in steady, quasi-one-dimensional flow. 



dpj 

dx 


f jA ; 


(306) 


Energy Conservation . - Now consider the energy balance for the same test 
mass of the jth gas discussed in the previous paragraph. Let ej be the internal 
energy per unit mass of the j“ gas and 4j be the total heat energy added to this 
gas per unit volume per unit time. qj is thus assumed to include all sources of 
energy addition to the j th gas in the test volume except for the work energy 
supplied by the partial pressure pj on the ends of the volume and by the body 
force fj . Thus , the energy balance for the test volume may be written as 


.x+Ax 

d / 1 


if 


+f - u j 2 ) A e Pj dx 


J qj a; dx + (Pj a; Uj) x - (Pj a; Uj ) x + Ax + f 


: + Ax 


fj Ag u j dx (307) 


Taking the limit of equation (307) as Ax— ►O and rearranging the result slightly 
yields the energy conservation equation for the j th gas in the form 



= (4j + fj Uj) A' e 


(308) 


where the specific enthalpy hj of the 


j* gas is defined by 


h i “ e i + 


(309) 
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Equations (302) , (306) , and (308) are the basic conservation equations which 
must be satisfied by the individual components of a gas mixture in steady quasi- 
one-dimensional flow. In general, the flow velocities uj of the individual species 
in these equations may be different; however, in the present analysis it is 
assumed that the flow velocity 

Uj - u 

(310) 

is the same for all species. For the heavy particles this should be a good 
approximation for pressures at least down to the order of 10"^ atmospheres, while 
for the electrons the condition (310) will be met provided there is no net flow 
of electrical current in the gas. It is further assumed that 

n n 


h S 'i -0 

(311) 

i-i i=i 


so that there is no net production of mass and no net body force on the flow. 

Then equations (302) , (306) , and (308) may be summed for all species j in the gas 
to obtain the usual quasi-one-dimensional conservation equations for the total 
flow 

— (pu A') - 0 
dx 

(312a) 

du dp 

^ dx dx 

(312b) 

d 1 

P u — (h +— u^) = q 
dx 2 

(312c) 

where 


i=i 

(313a) 

is the total density of the gas mixture. 


III 

aM- 

JV 

(313b) 
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is the total pressure 



(313c) 


is the total enthalpy of the gas mixture per unit mass, and 



(313d) 


is the total energy addition to the flow per unit volume per unit time from ex- 
ternal sources. For the present analysis, the only external source term to be 
considered will be the radiative energy loss from the gas, so that the source 
term (313d) becomes simply 

q - - q r (314) 

where q r is the radiated power per unit volume lost from the gas. 


For the most general type of one-dimensional flow, the species densities 
P j , flow velocities uj, and internal energies ej would all be independent vari- 
ables, so that the species mass, momentum, and energy conservation equations (302), 
(306) , and (308) would all be required to obtain the complete flow solution. 

Because of the assumption (310) made in the present analysis, however, these 
equations are no longer all independent, so that the body force |fj may be elimi- 
nated between the momentum and energy equations (306) and (308) to obtain the 
single equation 

£ (p i U A 'e + P i S' (U A e> - J " 2 Z (P > " A ' e> = A ' e <315) 

for the species internal energy ej . Further, it will be assumed in the present 
analysis that the exchange of translational energy among the heavy particles in 
the gas is sufficiently rapid to maintain thermal equilibrium among them, so that 
equation (315) is required only for the electrons. Thus, for a gas mixture con- 
taining n species, there are n + 3 independent flow equations (namely, the n species 
conservation equations (302) , the overall momentum and energy conservation equa- 
tions (312b) and (312c) for the mixture, and the electron energy equation (315) to 
determine the n species densities pj , the flow velocity u, and the electron and 
heavy-particle temperatures T e and T for the mixture. 


For use in the NATA code, it is convenient to express the above flow equations 


in terms of the species concentrations 

JL 

p w i 


in moles/gm, defined by 


r, - 


(316) 
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where Wj is the molecular weight of the j* species and p is the total density 
of the mixture defined by equation (313a). When this notation is introduced 
into equation (302) and the result simplified by the use of equations (310) and 
(312a), one obtains the species conservation equations in the form 


P u 


dyj 

dx 


W i 


Similarly, the electron energy equation (315) becomes 


P u 


— (7 e e e> 
dx 


Pe d P 


W p 


2 dx 


d y e 


dx 


= Se 


(317) 


(318) 


or, introducing the ideal gas equation of state for the electrons and noting that 
the flow energy term u 2 /2 will be negligible compared to the random thermal energy 
<r e for electrons because of their small mass. 


P u 


d / 3 \ d in p 

y e Rq T e ) - le R 0 T e 


dx \ 2 


dx 


(319) 


where Rq = 1.9872 cal/mole-°K is the gas constant per mole. 


To complete the specification of the problem, it is now necessary to specify 
the source terms rj , q r , and q e occurring in equations (312), (317), and (319). 

In general, these terms may be written as sums of contributions from the individual 
collisional and radiative processes occurring in the gas. Thus if the i* reac- 
tion process occurring in the gas is represented by the chemical formula (269) , 
then the total number of reactions of the type (269) occurring per unit volume of 
the gas per unit time in the forward and reverse directions will be respectively 


N fi = N o k fi ]Q (pyj) 1 " 1 ’ 
j = l 


and 


N ri = N 0 k ri I"I (Pr,> 11 

j=l 


(320a) 


(320b) 


where N 0 = 6.0225 x 10 ^^ ± s Avogadro’s number and k fi and k ri are respectively the 

forward and reverse reaction rates for process i in units of (mole/cm^)^ 1 sec 

The source terms r- s q r , and q e are then given in terms of the numbers of reactions 
(320) by i 



(321a) 



(321b) 


r 



i = 1 



N 0 




N r 


N n 


(321c) 


where the sum is over all reactions i occurring in the gas, /3ij = v ij - "ij is the 
net number of molecules of species j which are formed in a single reaction of type 
i in the forward direction, qfj and -qri are the mean energies lost from the gas 
by radiation in Nq reactions of type i in the forward and reverse directions, respec- 
tively, and and -e r j are similarly the mean energies gained by the electron 
gas in No reactions of type i in the forward and reverse directions. In general, 
the reaction rates k f j and k r j in equations (320) and the parameters qfj , q r j , (f- 
and e r j in equations (321) depend on the reaction under consideration and must be 
evaluated individually for each gas. The rate constants and parameters used in 
the standard gas models for argon and helium will be documented in Volume II of 
this final report (The NATA Code - User's Manual). 


Substitution of (321a) into (317) gives the species production equation (286) 
used in the chemical non-equilibrium model. Also, (312a) and (312b) are identical 
with equations (270) and (271), respectively. The set of governing equations for 
the electronic nonequilibrium model differs from that for the previously considered 
chemical non-equilibrium model in the following respects: 

(1) The electronic non-equilibrium model contains one additional dependent 
variable, the electron temperature, Tg. 

(2) There is one additional governing equation, the electronic energy 
equation (319) . 

(3) The global energy equation (312c) contains the radiative loss term (314) 
on the right, whereas the corresponding equation (272) in the chemical 
non-equilibrium model has zero in the right-hand side. Thus, in the 
electronic non-equilibrium model the total enthalpy decreases in the 
downstream direction. 


(4) The equation of state (273) is replaced by (313b), in which the partial 
pressure p. for the j th species is 


pj = R 0 Nj Tj 

Here Nj = pyj is the number of moles of 
and Tj is the translational temperature 
assumed to be equal to T for all of the 
and (322) give 


(322) 

the j th species per unit volume, 
for the species. Since Tj is 
heavy species, equations (313b) 


p = R 0 P (>e T e + yh T> 


( 323 ) 
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in which 



( 324 ) 


is the total molar concentration of the heavy species.* 


For use by the code, the species production equations are rewritten in the 
form (287) . The only modification of these equations in the electronic non- 
equilibrium model is that the equilibrium constant Kj in the expression (289) for 
the quantities x; ma Y t> e calculated from the electron temperature T e instead of 
the gas temperature T, for some of the reactions. In such cases, Kj is computed 
from 


K i = 


(R 0 T e ) 


ft 


exp 



(325) 


in place of equation (278) , 


Elimination of the derivative du/dx of the flow velocity between the con- 
tinuity and energy equations (312a) and (312c) gives 


dltx p 


dx 


d£n a; 


dx 


dh 


dx 


P u 


0 


(326) 


when (314) is taken into account. For the electrons, equation (292) is replaced 
by 

d H 


dx 


dT e 
'P e dx 


(327) 


If (291) is used to eliminate dh/dx from (326), and (292) and (327) are used for 
the heavy species and the electrons, respectively, there results 


'Ye ft 


dT e 

n 

n 


dT 

+ 

dx 

yiCpi 

dx JLvJ 


j=2 

) = 1 

9r 

d En0 

dfn Al 
2 e 


i 

H j dx 


(328) 


+ — - u 

pu 


dx dx 

This is the electronic non-equilibrium analog of equation (293) 
•Equation (324) is based on the NAT A convention that the electrons, if present, are species number 1 . 
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Differentiation of the equation of state in the form (323) gives 


dp dp / dT e dT \ 

— = R o- (y e T e + y h T) + Rqp. y y e >h -£) 


( dy e 

+ R ° P \ Te IT + T 



( 329 ) 


Substitution of (312c) and (329) into the momentum equation (312b) gives 


dh 

dx 


pu 


+ R 0 (ye + yh 1) 


+ Ri 


0 


d y e 


dyh 


dx dx 


d^n p 

dx 

)-° 


+ Rg 


dT„ 


dT 


17 + ^ 17 


(330) 


Elimination of dh/dx using (291), (292), and (327), followed by addition of (328), 
gives 


dT. 


dT 


dy, 


r ' ~ * »> * * T * 


dx 


e d Yh 

+ T 


Ye + Yh ^ 


Rr 


dx 


dx 

.2 


d(n Ag 
dx 


= 0 


(331) 


which is the electronic non-equilibrium analog of equation (299) . The electron 
energy equation (319) may be rewritten as 


2 dx +y * 

The total enthalpy, h Q 


d Te 

dx 

= h + 


Ye Rq T e 


d tnp 


= 0 


dx pu 

u 2 /2 , obeys equations (312c) and (314), 


i.e. , 


(332) 


dh 0 q r (333) 

dx pu 

To follow the changes in this quantity accurately, NATA treats hg as an additional 
dependent variable in the numerical integration. Thus, there are n + 3 dependent 
variables, T, T e , h 0 , and the yj , in the electronic non-equilibrium model. At 
each point x in the flow, the concentration derivatives dy;/dx are computed from 
equations (287), and the derivatives dT/dx, dT e /dx , and dh 0 /dx are obtained by 
simultaneous solution of equations (328), (331), (332), and (333).* The con- 
ditions at the point x + Ax are then calculated from the flow variables and their 
first-order derivatives at x, using the modified Runge-Kutta technique described 
in Section 7.5. Once T, T e , Hg , and they: have thus been determined at the new 
point, the specific enthalpy is computed from (275), which here takes the form 


h 


y e H e (T e ) + 



Yj Hj (T) 


(334) 


"Either d £np/dx or d £n A^/dx is also obtained from the simultaneous solution, but these quantities are not integrated numerically. 
Section 7.4 discusses the reasons for this procedure. 
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The flow velocity u Is then calculated from 


u = V 2 (h 0 - h) (335) 

Either the area A e or the density p is assumed known as a function of x . The 
other quantity is calculated from the continuity equation (243) . The pressure 
is calculated from the equation of state (323), the entropy from (251), and the 
Mach number from (300). 


7.1.3 Conservation of Chemical Elements 

Each of the individual chemical reactions (269) must be "balanced", i.e., 
must conserve the chemical elements. If this is true, the entire reaction system 
conserves the elemental composition of the gas mixture, regardless of the rates 
at which the individual reactions proceed. However, truncation errors in the 
finite-difference solution of the chemical rate equations lead to small changes 
in the elemental composition in each integration step. If allowed to accumulate, 
such errors could become objectionably large. In NATA, such gradual shifts in 
the elemental composition of the gas are prevented by adjusting the species con- 
centrations after each successful integration step. This adjustment is based on 
the following relations. 

The number of gram-atoms of the j th element per gram of the gas mixture is 
given by 


n 



i = 1 


-E 


(Ti) 


o x ) 


constant 


(336) 


In this equation, aij is the number of atoms of the j th element per molecule of 
the ith species, y j is the concentration of the i 1 * 1 species (moles/gm) at the 
current flow point, and (yj) 0 is the concentration in the upstream reservoir. 

The equation (336) can be solved to give the concentrations of the independent 
species in terms of the c- and the concentrations of the dependent species. (See 
Section 2.1.) From (336), 


c n 



i=l i=c+l 


(337) 


Multiplication of (337) by the inverse Ajjj of the square submatrix of «ij with 
i = 1 to c, followed by summation over * j , gives 
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ni- c ,k 


(338) 


c 



n 

Z 

i .5= r 4- 1 


for k = 1 to c, where ^i- c , k is defined by equation (6). 

In some earlier versions of NATA, the elemental composition was maintained 
by recomputing the concentrations y^ of the independent species, using (338), 
after every integration step. This procedure proved to be unreliable. If fhe 
concentration of one of the independent species is very much smaller than the 
concentrations of some dependent species containing the same chemical elements, 
the loss of accuracy on the subtractions indicated in (338) can become so severe 
that the independent species in question fluctuates wildly or even goes to a 
negative coiicentration. To avoid this problem, the present version of the code 
uses a more elaborate element-conservation algorithm which spreads the corrections 
over all of the species, more or less in proportion to their concentrations, in- 
stead of adjusting only the independent species. The adjustments Syj to the 
concentrations y< computed in an integration step are selected so as to minimize 
the sum of squares of the relative adjustments, 

»-EGt) ’ 

1 = 1 


subject to the constraint that equation (338) be satisfied by the adjusted con- 
centrations 

7i = y[ + (340) 

Prior to the adjustment, the concentrations of the elements have the incorrect 
values 


c i 


n 



i = 1 


(341) 


instead of the correct values cj based on the gas composition in the reservoir. 
Thus, before the adjustment 


c n 



j=l i=c+l 


(342) 
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Subtraction of (342) from (338) gives 




53 Sc )' A i k_ 2 Syi ‘ " i ‘ c - k 


i = 1 


i = c + 1 


where 

Sc ) a c i - c f 

Thus, the sum of squares (339) of the relative adjustments can be written 
terms of the adjustments for the dependent species: 


C ’8y x 2 


D(Sy c + 1 , . . . , Sy n ) 
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The conditions dD /d(8y- L ) = 0 for minimum D can be written in the form 
n 

y| a im d Ym = b i 
m = c + 1 


where 
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(344) 
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(346a) 


(346b) 
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In NATA, the system of linear equation (346) is solved for the adjustments to the 
concentrations Sy m of the dependent species using subroutine DSMS0L. The adjusted 
concentrations of the independent species are then computed using equation (338) . 
The entire element-conservation calculation is done in double precision to avoid 
excessive loss of accuracy in the final evaluation of the independent species 
concentrations using (338) . 


7.2 Method of Solution 

The basic method for calculating the non-equilibrium solution in NATA is 
numerical integration of the rate equations derived in Section 7.1. However, 
there are two difficulties which prevent immediate application of this method 
when the solution is started in the upstream reservoir: 

(1) Startup of integration near equilibrium . - The gas in the reservoir is 
assumed to be in equilibrium, with zero flow gradients. Also, for the 
kinds of nozzle profiles that are used in NATA, the derivative of the 
flow area, dA e '/dx, is finite everywhere; thus, the flow velocity is 
quite low in the far upstream region of the nozzle, and in consequence 
the gradients of the flow variables dT/dx , dp/dx , d y. /dx are quite small 
in that region, and the flow is close to equilibrium. For flows which 
are nearly in equilibrium, the numerical integration procedure is stable 
only for extremely small step sizes. Thus, it is impracticable to start 
the integration far upstream in the nozzle. 

(2) Flow upstream of the throat . - The gas conditions in the reservoir, to- 
gether with the nozzle geometry, determine what the mass flow must be. 
However, the correct mass flow for a given non-equilibrium flow problem 
is not known when the solution is started. The mass flow for the corres- 
ponding equilibrium flow problem is known, but typically differs slightly 
from the non-equilibrium value. If an attempt were made to compute a 
non-equilibrium solution by straightforward integration of the rate 
equations, assuming the equilibrium-flow value of the mass flow, most 
probably the solution would fail. Either the Mach number would reach 
unity before the throat was reached, in which case the solution would 
blow up because of the flow gradients becoming extremely large; or the 
Mach number at the throat would be less than unity, in which case the 
transition from subsonic to supersonic flow would be impossible. To 
obtain a valid solution in this way, an iteration to determine the 
correct mass flow would be required.* However, repeated numerical 
integrations of the flow equations from the reservoir to the throat 
would be extremely time consuming. 

These two problems are dealt with, in NATA, using methods which were evolved 
over a period of several years at Cornell Aeronautical Laboratory (refs. 1, 38). 
The difficulty in starting the integration near an equilibrium flow condition is 
avoided by treating the initial portion of the non-equilibrium solution by a 
perturbation method., in which the unperturbed solution is the inf inite-reaction- 
rate equilibrium flow. This perturbation method is explained in Section 7.3. 

The problem of calculating the non-equilibrium flow in the region upstream of the 
throat is handled using an inverse method which is documented in Section 7.4. 


•Methods of dealing with this problem are discussed by Bray (ref. 31). 
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Beyond the throat, and beyond the near-equilibrium region in which the perturbation 
method is employed, the solution is calculated by direct numerical integration of 
the rate equations, using the modified Runge-Kutta technique described in Section 
7.5. 

Figure 24 shows, in flowchart form, how these three methods are combined to 
generate the entire non-equilibrium solution. The solution is always begun using 
the perturbation method. In this method, the step from one flow point to the 
next is taken by decrementing the gas temperature, as in the equilibrium flow 
solution (Section 6.2). After each step, certain quantities SXj, which measure 
the departures of the reactions from equilibrium, are tested. When any one of 
these quantities reaches a specified size, the perturbation method is abandoned 
and the numerical integration is started. During the integration, the step from 
one flow point to the next is taken by incrementing the position coordinate, x . 

If the solution has not yet reached the "downstream region" (bounded by a point 
somewhat beyond the throat, at which the flow is already supersonic), the inverse 
method is used. Once the downstream region has been reached, the solution is 
generated using the direct method. If the downstream region is reached in the 
perturbation solution, then the inverse method is never used at all. 


7.3 Perturbation Method 

The flow in the far upstream portion of the nozzle, where the conditions 
differ only slightly from those in the reservoir, is nearly in equilibrium. The 
species concentrations and flow variables can therefore be expressed in the form 


Yj ~ Yj + $Yj 
T = T + ST 


(347a) 

(347b) 


P = P + S P 


(347c) 


in which the barred quantities are values for the equilibrium flow and 5yj , <5T , 
Sp are small perturbations due to departures from equilibrium. 


7.3.1 Perturbation of the Rate Equations 

A system of equations for determining the perturbations is needed. Such a 
set of relations can be derived from the chemical rate equations (287) and the 
flow equations, as follows. Substitution of (347) into the rate equations (287) 
gives, to first order, 

d y; d( Sy -. ) _ _ 

— + Pij (P i*i + X i Sp i + p i 5 *i> (348) 

i = 1 


The bar symbol in X; , P; denotes that a quantity is to be evaluated using the 
equilibrium flow values of the dependent variables, y- , T , and p . Thus, for 
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example, P- does not denote the value of Pj for the infinite-rate equilibrium 
flow, but rather Pi evaluated from (288) using the actual value of the rate 
constant kfi and the equilibrium values y- , T , and p . Equations (277) , (279) , 
(283) , and (289) thdn show that 5?i « 0 , so that_Xi is a measure of the departure 
of the i* reaction from equilibrium. Setting X; = 0 in (348) gives 


dyj d(Syj) 

dx dx 


2 sx i 

i =1 


(349) 



(350) 


Here, ( / dT)is the derivative dy./dTof (289), evaluated at equilibrium. 

From (289) , 1 
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P 



n ^ ik 

k =i 



(351) 


At equilibrium, = 0 or 
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K| (T) 
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k = 1 


y k 


ftk = ! 


Thus, (351) implies that 


d Xi 
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dK; 


\ dT f Ki(T) \ dT j 
From (278) , 

En Ki = — Pi (EnR 0 + En T) 
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dT 


[En Ki (T)] 


^ X) Pi ' Pi ° 
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(352) 
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Hence 


dfclK i ft y\ d /f*j° \ 

dT "t" 2m*i ' dT* \Rof / 

j =1 

Now, from equation (25) , 



Therefore, from (280), (353), (355), and (356), 



(355) 


(356) 


(357) 


Similarly, differentiation of (289) with respect to p and y. , followed by use of 
(352), shows that 1 



(358) 



Thus, (350) becomes 



(359) 


(360) 


From (349) and (360) , the entire right-hand side of (349) is of first order in 
the perturbations. Thus, dyj / dx is small of first order, and d ( <5yj ) / dx is there- 
fore of second order and can be neglected to within the accuracy of the pertur- 
bation calculation. Equation (349) is therefore rewritten as 



( 361 ) 
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Substitution of (360) then gives 



7.3.2 Equilibrium Derivatives of the Flow Variables 

The right-hand side of equation (362), dy./dx, denotes the gradient of the 
concentration of the j th species in the equilibrium flow. The equilibrium con- 
centration gradients depend upon dT/dx and dp/dx . Thus, all of the derivatives 
dyj/ dx for i = 1 to n, dT/dx , and dp/dx have to be determined. The relations 
(293) and (299), which are based on the one-dimensional steady flow equations, 
are valid for the equilibrium flow as well as for the non-equilibrium case. How- 
ever, the rate equations (287) cannot be used for the equilibrium case because Pj 
is infinite and x; is zero in this case. Thus, for the determination of the 
derivatives of the flow variables, n relations are needed to replace the n rate 
equations (287). The relations used for this purpose, in NATA, are the c element- 
conservation relations (336) , which give 


n 

d yi 

a ij IT = 0 (j = 1 ’-’ c) 

i =i 



(363) 


and the n-c equations (224) for the equilibrium mole fractions Xj of the dependent 
species. From equations (1) and (273), equation (224) can be written in terms of 
the equilibrium molar concentrations y. in place of the mole fractions Xj , in the 
form 1 ’ 


* w—mm 

fj = Kj (pR 0 T) V) ’" c_1 J"| yg V ' _c ’ e (364) 

1 = 1 


where Kj is the equilibrium constant (225) for the reaction forming the j th de- 
pendent species from the independent species. Logarithmic differentiation of 
(364) gives 
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(225) and 

(356), 
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( 365 ) 


(366) 


With the aid of equations (15) and (366), equation (365) can be rewritten in the 
form 


c 



Equations (293), (299), (363), and (367) are n+2 linear relations among the n + 2 
derivatives dj/j/dx , dT/dx , and dp/dx for the equilibrium flow. In NATA, this 
system of equations is solved using subroutine DSMS0L. 


7.3.3 Unperturbed Equilibrium Solution 

Equations (293) and (299) contain the derivative d£«A e /dx of the effective 
nozzle cross sectional area. During the non-equilibrium solution by the pertur- 
bation method, this geometric quantity is determined as follows. As in the NATA 
equilibrium flow solution (Section 6.2), the temperature is the independent vari- 
able. Successive points in the solution are generated by decrementing the 
equilibrium temperature, T . At each T, all of the equilibrium flow variables are 
computed just as in the equilibrium solution, using subroutine NEWRAP. From the 
equilibrium-flow values p and u of the density and flow velocity, the flow area 
ratio A e is calculated using the continuity equation in the form (265) . Then the 
geometric position x of the flow point in the nozzle is determined by solving 
equation (266) using subroutine FINDX. Finally, d fn A' / dx is calculated from the 
known nozzle geometry and the position x by calling subroutine GE0MAR.* 


*Note that din A' /dx = d In A /dx. 

e c 
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7.3.4 Equations Determining the Perturbations 

There are n+2 perturbations Syj , St, Sp to be determined at each flow point 
in the perturbation solution. Thus, n+2 independent linear algebraic relations 
among these perturbations are needed. Perturbation of the element conservation 
equations (336) gives c relations, 

n 

a ij = 0 

i =1 


Equations (362) are an additional n relations involving the Syj , but these re- 
lations are not all independent. Equations (362) were derived by perturbing the 
rate equations (287). As pointed out in Section 7.1.3, these rate equations 
automatically conserve the chemical elements provided the individual reactions 
(269) are balanced. Thus, the system of n + c equations 1362) and (368) contains 
the same information as the n equations (362) alone. Since the c equations (.368) 
are used*, only n-c of the equations (362) provide independent relations among 
the Syj . The equations (362) with j = c + 1 to n are the ones used. 

Another requirement for solvability of the equations for the Syj is that the 
reaction system provide independent chemical pathways for the formation and de- 
struction of all the species. If there is no reaction for forming one of the 
species, then its concentration cannot change and no choice of the Syj will allow 
the left hand side of (362) to match the equilibrium concentration gradient dy:/dx 
on the right. "If there are reactions involving a particular pair of species, but 
these reactions do not allow the two species to be formed or destroyed inde- 
pendently, then the species concentrations change in a fixed relationship and it 
is not possible for the left hand sides of equations (362) to match the independent 
equilibrium gradients dyj/dx 0 n the right. This requirement on the reaction system 
can be expressed in quantitative terms by considering the chemical reaction formula 
(269) in the form 


(j = 1. — . c) 


(368) 



(369) 


which may be obtained by combining (279) with (269) . In (369) , Mj denotes the j c 
species and j8jj is the number of molecules of this species produced in the i th 
reaction. Since each individual reaction conserves the chemical elements, t,he 
/3 r satisfy the c conditions for each reaction 



(k = l, , c) 
(i =1, .... r) 


(370) 


*An alternative approach would be to omit equations (368) and use all of the equations (362). The method used in NATA has the ad- 
vantage that equations (368) involve less computation and ensure accurate element conservation in the perturbation solution. 
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Equation (370) states that the net number of atoms of the k th element produced 
in the i* reaction is zero. Because of (370), only n-c of the ftj values de- 
fining a given ( i* ) reaction can be chosen independently; once these have been 
selected, the remaining values are all determined by conservation of elements. 
For example, if the ftj for the "dependent species" j = c + 1 to n (Section 2.1) 
are specified, then the ftj for the "independent species" j = 1 to c are all de- 
termined by element conservation. 


Now, if the species Mj are considered to define the directions in an n - 
dimensional vector space, then the ftj for each i are the components of a vector. 
The number of linearly independent reactions is then equal to the rank of the ft j 
matrix (ref. 39). According to (370), the rank of ftj can be no larger than 
n-c, regardless of how many reactions are included in the gas model. The rank 
can be smaller than n-c if there are too few linearly independent reactions, but 
in this case the perturbation solution will not work for reasons explained above. 
Thus, the reaction^ system is required to contain exactly n-c linearly independent 
reactions. In NATA, this requirement is applied by computing the rank of ftj 
using a standard technique (ref. 40). If the rank is found to be less than n-c , 
the case is terminated. Of course, all of the compiled-in gas models satisfy 
this requirement on the rank of ftj . 


Equations (368) together with equations (364) for j = c + 1 to n provide n 
relations among the n+2 perturbations Syj , St, dp. Two additional independent 
relations are required. One such relation may be obtained by perturbing the 
energy equation in its integral form (245) : 


8h + u 8u = 0 
From (275), 


(371) 


Sh = ^ t Hj ( T ) §7j + yjCpj ST] 
j =1 


(372) 


Perturbation of the continuity equation (243) gives 


p S u + u dp = 0 


(373) 


since the reactions do not affect the nozzle geometry. Use of (372) and (373) to 
eliminate Sh and Su from (371) gives 




( 374 ) 
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7.3.5 The Condition 8s = 0 

Equations (362) for j * c +1 to n, equations (368) for j= 1 to c, and 
equation (374) are a system of n+ 1 relations among the n + 2 perturbations S yj , 

ST, Sp . One additional relation is required to give a determinate system. 
Logically, this relation should be based on the momentum equation (244) or (271), 
which has not yet been applied to the determination of the perturbations. How- 
ever, the momentum equation is a differential equation; substitution of (347) 
into it would give a relation between the derivatives of the perturbations, 
whereas what is desired is an additional linear algebraic relation between Syj , 

ST, and Sp. To obtain such an algebraic relation, Cornell Aeronautical Laboratory 
(ref. 1) chose to introduce the approximation Sp~ 0 in the computer program from 
which NATA has been developed. In other versions of the same program, Cornell 
Aeronautical Laboratory used the alternative approximation that the perturbation 
in entropy is negligible, 

Ss = 0 (375) 

which is slightly more accurate (ref. 1). The current version of NATA uses (375) 
as the required additional relation between the perturbations Syj , ST, and Sp. 
Since both the infinite-reaction-rate equilibrium flow and the zero-reaction-rate 
frozen flow are isentropic (Section 6.2), it is plausible that the finite- 
reaction-rate non-equilibrium flow should be nearly isentropic in the initial 
region where the departures from equilibrium are still small. Thus, the pertur- 
bation in entropy from the equilibrium flow value, Ss , should be small in the 
region to which the perturbation technique is applicable. This supposition is 
amply confirmed by the results of integration of the rate equations using NATA. 

In general, the increase in entropy in non-equilibrium flow solutions is quite 
small, not only in the initial region but even in the region beyond the throat. 


The condition (375) can be expressed as a relation between ST, Sp, and the 
Syj as follows. From equations (1) and (296), the specific entropy (251) of the 
gas mixture can be written as 



R 0 En (p R 0 T) - Rq In Yj ] 


(376) 


To first order, the perturbation Ss is given by 
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(377) 
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From (376) 



Sj° ( T ) - Rfl [1 + In ( Rg p T yj ) ] 






[Cpj(T) - Rq] 



The final form of (378b) is obtained using the relation 

ii - 

dT T 


(378a) 


(378b) 


(378c) 


(379) 


which can easily be derived from equations (29), (31), and (356), and also fol- 
lows directly from the thermodynamic definitions of entropy and specific heat. 

The final form of equation (378c) is obtained using (296). Substitution of (378) 
into (375) and (377) gives 


u 

£ 

j = i 


Ss = / . * S i (T) ~ R 0 [1 + (RopT yj)]} S yj 




j [ s- r o ] m t 

j=i 


^o_ 
p w 


dp = o 


(380) 


7.3.6 Conditions for Starting the Numerical Integration 

Equations (368) for j =1 to c, equations (362) for j = c + 1 to n, 
equation (374), and equation (380) are the n+2 equations which are solved 
simultaneously for the perturbations Syj , ST , dp. These equations are set up 
in subroutine PERT and solved by calling subroutine DSMS0L. Once these pertur- 
bations have thus been determined, 5 X] is computed from equation (360) for each 
of the reactions. These S^ values are used to determine the point in the flow 
solution at which the perturbation technique is abandoned and the numerical 
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integration of the chemical rate equations is started. The rule for selecting 
this point is based on the following considerations: 

1) The accuracy of the perturbation calculation decreases with increasing 
distance down the nozzle, as the perturbations Syj , ST , dp become larger. 
The numerical integration must be started before the errors become ex- 
cessively large. 

2) Each .is an approximate measure of the departure of the corresponding 

reaction from equilibrium, since xj ~ X; + ^Xi = . When any one of 

the Sxi is very small,' the numerical integration is stable only for ex- 
tremely small step size, Ax . Thus, it is desirable to carry out the 
switchover to numerical integration as far downstream as possible, to 
avoid having to compute a large number of very small integration steps. 

The criterion for switching from perturbation to integration is 


'X 


S Xi 


12 C, 


(381) 


in which | d I is the largest absolute magnitude |8x; | for any of the reactions, 
max 

and C v is an input criterion value which is preset to 0.1. If | I turns out 
* 1 max 

to be greater than 1.2 Cy, at a point in the perturbation solution, a new point is 
computed for a temperature half way between this point and the preceding one. If 
necessary, the temperature interval is thus subdivided repeatedly in order to de- 
termine a point at which | S xi l max satisfies both of the inequalities in (381). 


Once a point satisfying the criterion (381) has been found, the perturbations 
8yj > 8T , 8 p are added to the corresponding equilibrium-flow quantities in accor- 
dance with equations (347), and the numerical integration is begun. 


Figure 25 illustrates the effect of the switching criterion on the non- 
equilibrium solution. Three NATA cases were run for the flow of air through 
standard nozzle 1, starting from reservoir conditions of 7000° R and 1 atm. In 
one run, C-^ had its standard value of 0.1; in the other two cases, values half 
and twice as large were used. Figure 25 shows the air temperature as a function 
of position along the nozzle, in the upstream region covered by the perturbation 
solution and the beginning of the numerical integration. The continuous curve 
shown in the figure is drawn through points of the perturbation solution. Each 
trace of unconnected graph symbols represents the results of the numerical inte- 
gration of the rate equations for one of the cases. The equilibrium and frozen 
flow solutions are also shown for comparison. The figure shows that increasing 
Cy, causes the code to begin the integration farther downstream. The differences 
between the three solutions provide an indication of the errors resulting from 
use of the perturbation technique. The vertical scale has been greatly expanded 
in figure 25 to show these small differences clearly. The solution for = 0.05 
is the most accurate, and should lie very close to the correct curve over most 
of the region shown. At the beginning of the integration for the =0.1 case, 
the temperature as given by the perturbation solution is too low by about 15 de- 
grees or 0.2 percent. For = 0.2, the corresponding error is 38 degrees or 
0.5 percent. After the integration has begun, each solution tends to relax 
toward the correct solution. For the cases shown in figure 25, at x= +1.27 cm 
the three solutions are separated by differences of about 1 degree out of 3770° K, 
or 0.03 percent. 
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7.3.7 Prevention of Premature Switch to Numerical Integration 

Por some reaction systems and some reservoir conditions, the [ <5y ; | values 
for different reactions can differ by many orders of magnitude. When this occurs, 
the criterion (381) would force the code to start the numerical integration at a 
point in the flow where some of the ISyjj values were extremely small. The re- 
sulting integration step size Ax required for stability of the finite difference 
equations would then be extremely small. The present version of NATA contains a 
computational artifice which circumvents this difficulty in most of the cases in 
which it would occur. At each step in the perturbation solution, the ratio 
R = I 8X: I . / I § X; I of the minimum and maximum I Sy. I values is computed. If R 

is smaller than an input DCHRAT, preset to 10 - 4, the code applies the following 
tests to the reaction i max giving the maximum ISXjl : 

1. Does any species j appearing in the reaction have a concentration yj 
which is less than or equal to a value GAMIN (preset to 10 - ^® mole/gm)? 

2. Is the reaction running in such a direction as to further reduce the 

concentration of this species; i.e., is 8: • P: xi <01 

max I max max 

If both of these conditions are satisfied, then the computational artifice is 
applied. It consists simply of increasing the rate constant for the i max th re- 
action by a factor of 1.1 x 10 ~Vr for the duration of the perturbation solution. 

This increase should make ISy. | equal to about 10^ times ISy: | , , which is a 

1 max min 

ratio allowing normal operation of the code. The justification for this artifice 
is as follows: First, it cannot affect the overall flow solution significantly, 

because it simply causes more rapid destruction of a species whose concentration 
is already negligible. Second, it is efficacious in most cases because the most 
common cause for an exceptionally high | Sy; | value is the presence of a species 
of high formation enthalpy Hg° whose equilibrium concentration is very low because 
the temperature is relatively low. Under these circumstances the equilibrium con- 
stant (225) for formation of the species, which contains a factor e -Ho°/ R oT £s 
very sensitive to the temperature because the exponent Hq°/r 0 t is large. Thus, 
the equilibrium flow solution demands a rapid decrease in the species concentra- 
tion. The inability of the actual reaction rate to follow this rapid equilibrium 
change is responsible for a large departure of the reaction from equilibrium and 
thus a large | Syj | . For example, this situation arises when a gas model including 
Ar, Ar + , and the 3-body recombination reaction 

Ar + + e — + e~ -> Ar + e~ 

is used at relatively low temperatures where the equilibrium concentration of 
Ar + is very low and the electron concentration may also be low. 


7.3.8 Neglect of Electronic Non-equilibrium Effects 

In runs based on gas models including electronic non-equilibrium, the in- 
equality of the electron and gas temperatures and the loss of energy by radiation 
are neglected in the perturbation solution. These features are "switched on" 
together with the numerical integration. The neglect of temperature non-equilibrium 
in the perturbation solution appears to be a reasonable approximation, since the 
gas is assumed to be nearly in equilibrium in the region where the perturbation 
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solution is used. The neglect of radiative losses is a poorer approximation, 
since the radiated power per unit volume is greatest in the high-temperature, 
high-density region near the reservoir. However, because of geometric approxi- 
mations used in NATA, inclusion of the radiative losses, in the perturbation 
solution would give radically incorrect results. In the upstream region, for 
x -» - oo , NATA nozzle profile curvefits open out in conical fashion rather than 
going to infinite area at a finite distance upstream of' the throat. Thus, the 
reservoir is, in effect, assumed to be an infinite distance upstream. Since the 
flow velocity is finite everywhere, the transit time of an element of gas from 
the reservoir to the throat is also infinite. This geometric idealization causes 
no difficulties in adiabatic flows. However, if the gas were allowed to emit 
radiation for an infinite time, it would lose all of its energy. Thus, the 
neglect of radiative losses in the perturbation solution allows the code to 
simulate the actual flow of a radiating gas which is heated in a region a finite 
distance upstream of the throat. 


7.4 Inverse Method for the Upstream Region 


The perturbation method used to start the non-equilibrium solution is based 
upon the equilibrium solution, and thus assumes the equilibrium-flow value of 
the sonic mass flux, m* (Section 6.2). When the numerical integration of the 
rate equations is started, the initial density p, velocity u, and effective area 
ratio A e at the switchover point thus have values which are consistent with the 
equilibrium sonic mass flux: p u A e = m, . If the integration is started down- 
stream of the throat, in the supersonic portion of the flow, it is carried out 
in a straightforward manner as outlined in Section 7.1. However, if the inte- 
gration is started upstream of the throat, a straightforward numerical integra- 
tion of the flow equations would almost certainly fail in the throat region 
because the sonic mass flux required to allow a smooth passage from subsonic to 
supersonic flow at the throat would differ slightly from the equilibrium value. 
To avoid this difficulty, NATA uses an inverse procedure when the integration is 
started upstream of the throat. This inverse method assumes, on the basis of 
previous studies (ref. 41), that the non-equilibrium density distribution up- 
stream of the throat and the sonic mass flux differ only slightly from the 
equilibrium flow values. The p(x) and m* for equilibrium flow are assumed to 
determine -the flow, in place of the specified nozzle geometry. To obtain a 
smooth representation of p(x) with smooth derivatives, an analytical curvefit 
to p(x) is used in place of the data at discrete flow points provided by the 
NATA equilibrium solution. The form of the curvefit is based on the approxi- 
mation of isentropic flow of a perfect gas with constant specific-heat ratio, 
y, In any standard elementary text on aerodynamics (ref. 42), the following 
relations are shown to apply to such a flow: 



pu A e 


y 

p y p ° 

(382a) 

y - 1 

p y- 1 p 0 




(382b) 

P0 y 


(382c) 


= p„ u„ = constant 


- 147 - 



P* ( 2 

\ i/( r — i ) 

(382d) 

r — < 

+ 

II 

1 £ 

) 

u. = = 

o 

CL 

r*s 

(N 

(382e) 

* * 

v y + 1 p o 


where subscript 0 

denotes reservoir conditions and * denotes 

sonic conditions. 

Elimination of u ; 

p , and u* from these four equations gives a 

relation between 

the density and the effective area ratio: 
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The value of a Is determined from equations (382d) and (384a) , 
P*' “ 


(a + 2) 


PO 


which give 


(385) 


The ratio / p 0 is known from the equilibrium solution. Equation (385) is solved 
for a using the Newton-Raphson method. The object of the iteration is to make 

/ » d 
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From the n estimate, a n , of a , the (n+1) estimate is calculated as 
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The derivative in the denominator of (387) is obtained by evaluating 
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(388) 


at a = a n . The iteration is continued until both |g | and l« n + i-a,i! ate less than 

or equal to 10~5, Then C is evaluated by substituting a and p„/p 0 into (383). 

»» 

Figure 26 illustrates the accuracy with which equation (383) fits the re- 
sults of actual equilibrium-flow calculations. The flow problem is the same one 

on which figure 25 was based. The curve represents equation (383), and the points 

are results of NATA equilibrium calculations. 


The non-equilibrium solution by the inverse method proceeds as follows 
each new value of x in the numerical integration 


At 


the geometric area ratio A (x) 
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is calculated by calling subroutine GE0MAR. Even if the run is one in which the 
boundary layer on the nozzle wall is to be included, the effective area ratio A e 
is assumed equal to the geometric area ratio A during the solution by the in- 
verse method. Then the density ratio p/ pg is calculated by solving equation (383) 
numerically. This computation is performed in subroutine GE0M. Thus, p/ Pg is de- 
termined as a function of the position x in the nozzle. The derivative dl n p/dx 
is also calculated from the relation 


d £n p 

dx 



2 C 
2 

- (a + 2)C 


d En Ag 
dx 


(389) 


which can be derived by taking the logarithm of (383), differentiating, solving 
for the derivative, and using (383) to re-express the result in somewhat simpler 
form. 


The flow equations (293) and (299) for the chemical non-equilibrium model, 
which are integrated during the non-equilibrium solution, contain both d in p/ dx 
and d£nA e '/dx . The same is true of the equations (328) and (331) for the elec- 
tronic non-equilibrium model. Either of these derivatives could easily be elimi- 
nated from the equations. Both are retained in order to accommodate both the 
normal, direct integration and the inverse procedure in the same set of formulas. 
In the direct solution, dEnA' e /dx is obtained from the nozzle geometry and dEnp/dx 
is determined by solving the flow equations. In the inverse method, d£n p/ dx is 
determined from the given nozzle geometry and the area-density relation (383) , as 
described above, and d£nA e '/dx is determined by solving the equations. 


Neither dEnp/dx nor dEn A e ' /dx is actually integrated in the Runge-Kutta 
routine (KNKT, Section 7.5). The quantities obtained from the integration at 
each flow point are T and yj (and T e and h Q in the electronic non-equilibrium 
model). From these quantities the enthalpy h is computed using (275). The 
velocity u is obtained from the integral form of the energy equation (245) . The 
continuity equation (265) is then used to calculate the density p in the direct 
solution or the area ratio A e in the inverse procedure. The effective area 
ratio A e thus determined in the solution by the inverse method is slightly dif- 
ferent from the actual area ratio at the flow point, A g (x) , which was used in 
determining the density from (383). This difference results from deviations of 
the non-equilibrium temperature and concentrations from the corresponding 
equilibrium values, which produce a difference in the calculated flow velocity. 
Typically, A e is a little smaller than A g (x) and drops below unity in the vicinity 
of the throat. f 

A e has its minimum at the point where the non-equilibrium flow velocity is 

equal to a = \J{ dp /dx) / (dp /dx ) . This point is near, but not necessarily at, the 
geometric throat of the specified nozzle geometry. In a reacting gas, there is 
a second sound speed besides a, the so-called "frozen" sound speed, c. The point 
in the nozzle at which the flow velocity is equal to c is a branch point; beyond 
the branch point, the governing equations permit two solutions, a supersonic 
solution (which is the one desired) and a subsonic solution. In nozzle flow of 
a non-reacting gas, a and c are equal and the branch point and the point of mini- 
mum area both lie at the geometric throat. In a reacting gas, the frozen sound 
speed is a little larger than a (rdfs. 31, 41), and thus the branch point is 
slightly downstream of the point of minimum area. 
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In NATA, the solution by the^ inverse method is continued beyond the throat 
until the calculated area ratio A e reaches a value (1.01) which typically occurs 
downstream of the branch point. Then the specified nozzle geometry is adjusted 
to make the actual effective area ratio A e continuous with the value X e computed 
in the inverse-method solution, and the solution is continued by direct inte- 
gration of the rate equations. The adjustment of the nozzle geometry is carried 
out by setting an area rescaling factor, R a , equal to the value of 



at the switchover point from the inverse method to the direct solution. This 
factor typically has values in the range 0.95 to 1.00. 

After the direct integration has been begun, the value of d£«p/dx computed 
at each point is checked for sign. If dfn p/dx goes positive slightly downstream 
of the switch point, it is assumed that the program has started to follow the 
subsonic downstream solution instead of the desired supersonic solution. This 
condition could arise if the branch point were actually downstream of the switch 
point. To recover the correct solution, the program restarts the inverse-method 
solution at the switch point, using previously stored data, and extends it farther 
downstream to a new switch point. If the d£«p/dx > 0 condition is encountered 
well downstream of the switch point, or if more than four restarts prove to be 
required, the case is abandoned. 


7.5 Numerical Integration 


7.5.1 Treanor Integration Technique 

The governing equations (287), (293), (299), (328), (331), (332) for the 
non-equilibrium nozzle flow problem are first-order differential equations of 
the form 


dy 

dx 


f (x, y) 


(391) 


Here the dependent variable y represents a concentration y: , the gas temperature 
X, the electron temperature T e , or the stagnation enthalpy h 0 , and x is the axial 
position coordinate in the nozzle. Actually there are n+1 or n+3 dependent vari- 
ables y and functions £ , and each function f depends upon all of the y’s . The 
subscripts needed to distinguish the different variables and functions are omitted 
in (391) to simplify the notation in the following discussion. 

Numerous numerical integration methods for systems of equations like (391) 
are available (ref. 43). The earliest version of the Cornell Aeronautical 
Laboratory program from which NATA was derived used a fourth-order Runge-Kutta 
method based on the equations 


Ay 


1 _ 

6 


Ax ( f j + 2 (2 + 2 f 3 + ) 


(392) 
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where 


h = f ( x i> 7i> 



(393) 

i 

*2 = X 1 + 7 Ax 

y 2 = yi + 

1 

— ft Ax 
2 1 

(394a) 

1 

Xj = x, + — Ax 

3 1 2 

y 3 = n + 

1 

— (j Ax 
2 2 

(394b) 

X^ = Xj + Ax 

Y4 =71+ 

fj Ax 

(394c) 


Here is a value of x at which the -tfalue y^ of y is already known, and y\ + Ay 
is the computed value of y at x = x^ + Ax . This integration algorithm is 
generally satisfactory. The truncation error in one step is of order (Ax)^, 
so that high accuracy can be achieved by suitable choice of the step size. Ax . 
Like other explicit integration methods, it is conditionally stable, i.e., 
stability can be maintained by using a sufficiently small Ax. However, a 
special problem arises in the application of this fourth-order Runge-Kutta 
technique to the non-equilibrium flow problem. The Pj factors (288) for the 
various reactions can differ by several orders of magnitude because of dif- 
ferences in the rate constants k fi and the concentrations of the partici- 
pating species. Consequently, some of the species concentrations may relax 
toward equilibrium very rapidly while others do so very slowly. For a species 
which is near equilibrium, the rate equation (287) takes the approximate form 


dy 

dx 


P(y-y) 


(395) 


where y is the local equilibrium concentration and P is an inverse relaxation 
distance. Differentiation of (395) shows that each derivative of y is P times 
larger than the derivative of previous order. Thus if P is large (i.e., if the 
relaxation distance is small), the higher order derivatives d n y/dx“ are very 
large. Under these circumstances, the Runge-Kutta technique is unstable except 
for very small step sizes Ax such that P • Ax <0(1). 


To deal with this problem, NATA uses a modification of the fourth-order 
Runge-Kutta technique developed by Treanor (ref. 44). When one of the P’s is 
large, Treanor' s technique normally allows the use of a much larger step size 
than is possible in the Runge-Kutta method. When the p’s are all small, it 
reduces to the Runge-Kutta method. 

Treanor' s technique is based on the assumption that, over the interval Ax, 
equation (391) can be approximated in the form 


— = f(x, y) = - P (y-yi ) + A + B (x«- ) + — C (x- Xi ) 2 


(396) 
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This relation can be integrated analytically from x = xj to x = xj + Ax to give 
the change in y over the interval: 


Ay = Ax tA Fj + (B Ax) F2 + (C Ax 2 ) F 3 ] 


where 

F 1 

f 2 

F 3 


-PAx 

1 — e 
PAx 

-PAx 

P Ax — 1 + e 
(PAx) 2 

(PAx) 2 - 2(P Ax) + 2 - 2 e 
2(P Ax) 3 


(397) 


(398a) 

(398b) 

(398c) 


The four coefficients P,A , B , C in (396) are determined by fitting this relation 
to data at the following four (x, y) points: 


= X 1 

y = yi 

(399a) 

1 


(399b) 

= Xi + — Ax 
1 2 

y = y 2 


1 

= xj + — Ax 

y = y 3 

(399c) 

= xj + Ax 

y = y 4 

(399d) 


The yj will be chosen at a slightly later stage in the analysis. Substitution of 
(399) into (396) gives four linear equations for P , A » B , C . The solution of 
these equations gives 



y$ - y 2 


(400a) 


A = t x (400b) 

BAx = — 3 ( f x + Py* ) + 2 ( f 2 + Py 2 ) + 2 (f ? + Py 3 ) - (f 4 + Py 4 ) (400c) 

C Ax 2 =4 [(f x + Pyj.) - (f 2 + Py 2 ) - (f 3 + Py 3 ) + (f 4 + Py 4 )] (400d) 

Substitution of (400) into (397) gives 
Ay = Ax Ifj FJ 

+ [ - 3 (f x + Pyi ) +2(f 2 + Py 2 ) +2(f 3 + Py 3 ) - (f 4 + Py 4 ) 3-F 2 (401) 

+ 4 [(fj + Pyj) - (f 2 + Py 2 ) - (f 3 + Py 3 ) + ( f 4 + Py 4 )3 F 3 I 
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The quantities y 2 , f and y3 are chosen, to be the values given by (394a) and 
(394b). Based on these values, the formula (400a) for P becomes 


P = _ 


2 

17 


— fi 


(402) 


When(P Ax) is large, the y 4 value given by (394c) can be considerably different 
from the correct value on the solution curve. To avoid the errors in B and C 
which would result from the use of (394c), a better approximation to y 4 is ob- 
tained by evaluating (397) with C = 0; i.e.. 


y 4 - Yi + Ay = y 1 + Ax [ A' Fj + (B' Ax) F 2 1 


(403) 


where the coefficients A 'and B'are determined by fitting equation (396) to the 
three points (xj , y 4 ) , (x 4 + Ax/2 , y 2 ) , and (xj + Ax/2 , y 3 ) . The expression (402) 
for P remains valid, since it depends only on the evaluation of y 2 and y 3 . A 
and (B' Ax) are found to be 


A' = fj 

B' Ax . 2 [(f 3 + Py 3 ) - (fj + P yi )] 

Hence, (403) gives 

y 4 = y x + Ax [2 f 3 F 2 + f 1 (Fj - 2F 2 ) + f 2 (P Ax) F 2 ] 


(404a) 

(404b) 

(405) 


This result is used, in place of equation (394c), for evaluating (B Ax) and (CAx 2 ) 
from (400) and in the final integration formula (401). 


For (P Ax)-. 0 , equations (398) show that 


Eim 

PAx-* 0 
tim 

P Ax -> 0 
Elm 

P Ax -> 0 


F, = 1 


F 3 = T 


(406a) 

(406b) 

(406c) 


Hence, for (P Ax) -> 0 , equation (405) for y 4 reduces to (394c), and the Treanor 
integration formula (401) reduces to the fourth-order Runge-Kutta formula (392). 
Thus, in a region of the integration where (P Ax ) is small, Treanor’s technique 
behaves like the fourth-order Runge-Kutta method, which is known to perform well 
under such circumstances. On the other hand, for large (P Ax) , Treanor’ s tech- 
nique allows the use of a much larger step size without instability. 

Equation (402) specifies the evaluation of P as the ratio of two differences. 
If ^ , f 2 , and f 3 are nearly equal, the loss of accuracy on the indicated sub- 
tractions can produce large errors in P, or even a negative value of p. However, 
when these slopes are nearly equal the fourth-order Runge-Kutta method is satis- 
factory. Thus, if the magnitude of (f 2 - ^ ) / fj is found to be smaller than 10-4, 
or if P is computed to be negative, the Runge-Kutta formula is used. 
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The preceding discussion has been phrased in terms of the solution of a 
single differential equation. The lion-equilibrium flow problem, of course, in- 
volves a system of n+1 (or n+3) coupled equations for the concentrations yj and 
temperature T^(snd for T e and ho in the case of an electronic non-equilibrium 
model). In|^5§4^ Treanor's integration technique is applied to each of these 
differential equations individually. This procedure is the same as the one used 
when the Runge-Kutta method is applied to a system of differential equations. 

Lomax and Bailey (ref. 45) have studied the stability and accuracy of 
several integration methods, including Treanor's, when applied to the equations 
of steady, one-dimensional, non-equilibrium flow. Their principal findings with 
regard to Treanor’s method were as follows: 

1. When applied to asingle differential equation, the method is stable 
for large step sizes, provided the parameter P is sufficiently close to 
the eigenvalue of the linearized equation. 

2. When the method is applied to a system of coupled equations, its stability 
depends on the coefficients. It works very well if the equations are 
nearly uncoupled, but the coupling is difficult to assess a priori . The 
allowable step size for a system of coupled equations depends upon the 
eigenvalues of the system. 

3. With proper step size control, the method is stable and gives accurate 
results. 

On the whole, Lomax and Bailey preferred an implicit integration technique 
to Treanor's method. However, the advantages of the implicit technique were most 
apparent in non-equilibrium shock wave calculations, in which the system under 
consideration approaches an equilibrium state downstream. Treanor's explicit 
integration formula has to use a small step size in near-equilibrium regions. In 
the nozzle flow problem, the departures from equilibrium generally increase in the 
downstream direction, so that the step size can be increased as the calculation 
proceeds. Thus, Treanor's method is expected to be more satisfactory for nozzle 
flow-problems than for shock wave problems. 


7.5.2 Step-Size Controls 

Proper control of the step size is extremely important in any numerical in- 
tegration using an explicit finite difference scheme. If the step is just a 
little too large, the solution becomes noisy and inaccurate. If the step is in- 
creased still further, the solution becomes unstable and small errors are ampli- 
fied without limit. On the other hand, if the step is much smaller than necessary, 
the solution requires a needlessly large number of steps and wastes computer time. 
In nozzle flow calculations, the step must be small initially, where the gas mix- 
ture is still near an equilibrium state, but becomes stable for much larger steps 
far downstream where the reactions are nearly frozen. Thus, a procedure for de- 
termining and maintaining the proper step size in each part of the calculation is 
needed. 

The step size controls in NATA are based on those described in ref. 1 for 
the Cornell Aeronautical Laboratory program from which NATA was derived. The 
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controls in NATA are' somewhat more elaborate and restrictive to meet the needs 
of the electronic non-equilibrium models, which are distinctly less stable than 
conventional chemipal non-equilibrium models. The controls are based on tests 
for the acceptability and accuracy of the computed results, both at intermediate 
points in each Treanor-Runge-Kutta integration step and at the end of each step.. 
The tests at the end of each step in conventional chemical non-equilibrium models 
are as follows: 


yj_>0 (j=l, - 

n) 

(407a) 

d £n A e 

X • : > - 0.01 

dx 


(407b) 

Ayj 

| | < GTEST (j = l, .. 

Tjb ~ 

n) 

(407c) 

i AT , 

< TTEST 

T , _ i — 


(407d) 


Here the subscript b denotes values at the beginning of the integration step. 

The criterion values GTEST and TTEST are under input control. They are preset 
to the values GTEST = 0.1 and TTEST = 0.05. 

The first test, (407a), requires that the concentrations be positive. The 
second requires that dAg/dx be negative upstream of the throat and positive down- 
stream of the throat; this test is effective only during the solution by the in- 
verse method, in which p (x) is assumed given and d in A e / dx is computed from the 
flow equations. The small negative value (-0.01) on the right allows for pos- 
sible slight displacement of the non-equilibrium sonic point from the geometric 
throat. The last two tests (407c, d) set limits on the changes in the dependent 
variables in any one step. 


In the case of an electronic non-equilibrium model, the following additional 
tests are performed at the end of each integration step: 


|— >| 

| <. TETEST 


(408a) 

T eb 

Abo 

<_ HTEST 


(408b) 

Nib 1 




1 A^ | 

< max QTEST • 

1 ^eb 1 ’ ®qm ^ 

(408c) 

where TETEST, HTEST , and QTEST are inputs preset to 0.05, 

0.01 and 0.1, respec- 


tively, and Dq m is the maximum value of QTEST • | q e b | computed previously in the 
case. The tests (408a) and (408b) set limits to the changes in the electron 
temperature T e and the total enthalpy ho in each step. The test, (408c) controls 
the fluctuations in the energy transfer q e to the electron gas. This test is 
needed in the electronic non-equilibrium model because in the upstream region, 
where some of the reactions are near equilibrium, 4e is extremely sensitive to 
small changes in the temperatures T and T e and some of the mole fractions yj . 

In turn, q e affects the electron temperature through equation (332) . The 
sensitivity of 4g to the dependent flow variables is the principal cause of the 
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poor stability of the electronic non-equilibrium models, which requires the use 
of small step sizes in non-equilibrium flow calculations based on such models. 

The form of the test (408c) allows for the possibility that q e may pass through 
zero at some point in the solution. 

The tests performed at intermediate points in the Treanor-Runge-Kutta inte- 
gration step are generally more lenient than those done at the end of the step. 

In the chemical non-equilibrium models, the tests at intermediate points are (407a), 
(407b) , and 


T > 0 


(409a) 


H 'C hn 


(409b) 


In electronic non-equilibrium models, the following additional tests are per- 
formed at intermediate points. 


T e _> 0 


h 0 2 0 


A 4e I l. 2 D qm 


(409c) 

(409d) 

(409e) 


At the start of each integration step, the values of yj , T , T e , x, hg, A e , p, 
and s are saved in a separate set of storage locations. If any of the conditions 
(407) to (409) is violated at any point where it is applied, the flow variables 
are reset to their values at the start of the step, the step size is reduced, and 
the step is repeated. If the step fails again, this procedure is repeated until 
a successful step which passes all of the tests is obtained, or until the step 
size falls below lO - -^® cm. If the step size goes below this limit, or if a 
successful step has not been achieved after 30 tries, diagnostic data are printed 
out and the case is abandoned. 

Changes in the step size are based on two factors, SC (used for increasing 
the step) and SCD (used for decreasing the step). SC is initially 1.1. If NQS 
successful integration steps have been taken without the need for a step size re- 
duction, the step size Ax is multiplied by SC. If the step has thus increased 
NQS times without the need for a reduction, SC is increased by 0.1. NQS is an 
input (with input name NQSI) , preset to 4. Thus, so long as the integration 
proceeds with no failures of the tests, the step size is increased repeatedly, 
and the factor SC by which Ax is multiplied is also increased. 

When a failure of one of the tests occurs, the following actions are taken: 

1. The counters for successful steps and successful increases in step size 
are both reset to zero. 

2. SC is decreased by 0.1 (but not below 1.1). 

3. Ax is divided by SCD. 

4. SCD is multiplied by 1.1. 


5. The step is restarted from its beginning. 
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If repeated failures of the same step occur, the increase in SCD on each attempt 
allows the program to reduce the step size very sharply in a small number of 
cycles. However, once the step has been taken successfully, Ax is reset from its 
(possibly very small) value used in the integration to max(Ax , 0.7 Ax ol( j) # where 
A x old is the step size used in the previous successful step. Also, after a 
successful step, SCD is reset to the current SC value. 

The initial step size at the beginning of the numerical integration is an 
input (DELTXI). It is preset to 0.01 cm. However, if 

Ax m = 100 |aYii min ( 4 io) 

is less than the input value of Ax, the step size is reset to Ax m . This provision 
takes account of the smaller step size required when the switch from the pertur- 
bation technique to numerical integration occurs in a region where one or more of 
the reactions are still very close to equilibrium. 

The preset values for TTEST and GTEST do not imply that the program allows 
errors of 10 percent in the concentrations or of 5 percent in the temperature. 

In fact, the cumulative truncation errors' which occur when the step size is con- 
trolled by (407) are typically of the order of 0.1 percent or smaller. 


7.5.3 Freezing of Minor Species 

Previous versions of NATA were subject to a mode of failure which typically 
occurred far downstream of the nozzle throat. The symptoms of this type of 
failure were as follows: The concentration of a minor species, already very low, 

would begin to decrease more and more rapidly, until changes in the concentration 
of this species were controlling the step size in the non-equilibrium integration. 
The decrease in concentration would continue to accelerate, forcing successive 
decreases in the step size. Thus, the rate of progress of the solution would be- 
come quite low. In some cases, the step size would become so small that loss of 
accuracy on subtractions would lead to major fluctuations in some of the flow 
quantities, and the code would be unable to proceed with the non-equilibrium 
solution. 


The mechanism of this type of behavior is as follows. The gas models used 
in NATA contain numerous reactions of the type 

S + A B + C (411) 

Equation (283) for the rate of change of species S due to this reaction is 
dy s 

— = (kfyS/A - k rTB I'd (412) 

ax u 


since iq=v;'= 2 for the reaction (411). Now, the reverse rate constant k r is 
given in terms of kf and the equilibrium constant K by equation (277) . Hence 
(412) can be rewritten 


j_ £ys k fp / xb re \ 

y s dx u V A Ky s / 


(413) 
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The equilibrium constant (278) for the reaction (411) is 



Now, assume that the reaction (411) has been written so that products (B and C) 
are much more stable than the reactants (S and A) at low temperatures. Then the 
exponent in (414) is positive. Because of its factor 1/T, this exponent can- be- 
come large in magnitude as the temperature falls. In such a case, the equilibrium 
constant K itself can become exceedingly large. Moreover, under such circum- 
stances, K increases very rapidly as the temperature decreases. Even though y s 
is decreasing, the increase in K can be so rapid that the second term on the right 
in (413) becqmes and remains small in comparison with the first term. If the con- 
centration yA of the second reactant A is large compared with ys , the reaction has 
little effect on yA . Under these conditions, dPnyg/dx is negative and varies 
only slowly with x . Thus, yg approaches zero in an approximately exponential 
fashion, 

y s — e~ x ^ L , L = u/k fP y A (415) 

The type of failure described above would occur under these circumstances if the 
characteristic relaxation length L were small. For some of the reactions used in 
NATA gas models, the rate constant kf is sufficiently high to produce such failures 
in some cases. 

The current version of NATA contains a procedure designed to prevent failures 
of this type. Whenever the following three conditions are satisfied: 

(1) An excessive change in the concentration yj of a species, over a complete 
integration step, has forced a reduction in the step size; 

(2) The concentration yjt, of the species at the start of the integration step 
was less than an input value GAMIN, preset to 10 _ 10 mole/gm; 

(3) The concentration yj decreased during the step, 

then the concentration of the species j is frozen by switching off all of the re- 
actions which produce or destroy the species. The rationale for this procedure 
may be outlined as follows: 

(1) It prevents failures of the type described above by switching off the 
reaction causing the rapid destruction of the minor species S (along 
with the other reactions affecting the concentration of the species). 

(2) Because the concentration of S is very low, and since the net effect of 
the reaction system is to transform S into other species, freezing the 
concentration of S cannot significantly affect the concentrations of the 
other species during the remainder of the non-equilibrium flow solution. 

When this procedure is used, NATA prints out a message identifying the species 
whose concentration is being frozen and the reactions which are being switched off. 
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7.6 Boundary Layer Effects 


The method used in the approximate calculation of laminar boundary layer 
development on the' nozzle wall has been documented in Section 5 . The boundary 
layer calculation for the non-equilibrium solution follows an approach similar 
to that used for the frozen and equilibrium calculations (Section 6.3). In the 
upstream and throat regions, the boundary layer is computed but is not coupled 
to the inviscid flow. A short distance beyond the throat, the coupling is "switched 
on". The motivation and justification for this approach have been explained in 
Section 6.3. 

In the frozen and equilibrium flow solutions, the coupling of the inviscid 
flow to the boundary layer is turned on at the third point beyond the geometric 
throat. In the non-equilibrium solution, the uncoupled solution is carried a 
little farther downstream, to the point at which the switch from the inverse 
solution to the direct integration is made. This procedure allows the effect 
of the displacement thickness on the area ratio at the switch point to be included 
in the area rescaling factor R a equation (390). Also, it avoids stability 
problems which might result from switching the coupling on very close to the 
throat, where the solution is extremely sensitive to the area ratio. 


In the uncoupled region upstream of the switch point, the effect of the 
boundary layer upon the effective area ratio A e is neglected ; i. e. , A e is assumed 
equal' to the geometric area ratio, A g (x). Downstream of the switch point, A e 
is calculated from A g (x) and the displacement thickness 8 using equation (126), 

(130), or (134) depending upon the type of nozzle geometry. Also, d In A^/dx = d £n Ag/dx 
in equations (293) and (299) is calculated from the derivative of the appli- 
cable area ratio formula. The derivative dS/dx appearing in these formulas is 
evaluated by the simple first-order difference expression 



in which 8 and x are the displacement thickness and axial coordinate at the current 
flow point, and Sq » x o the corresponding values at the preceding flow point. To 
avoid an abrupt ^discontinuity in d In A e /dx a t the point where the coupling is 
switched on, d 8 /dx i s built up gradually from zero over 29 integration steps, 
according to the formula. 



where 


w 


1 

30 - i 


(418) 


and i is an integration step counter which is initialized to zero in the step in 
which the coupling is switched on. Also, (d8*/dx) 0 is the final dS*/dx calculated 
in the preceding step using (417). If (d8Vdx) c were constant, equations (417) 
and (418) would give a linear variation of dS*/dx with the counter i, from 0.0333 
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(dS*/dx) c at 1 = 0 to 0.9667 (dg*/ dx) c at ‘ = 28. For i > 28, equation (418) is 
no longer used. Instead, the weight factor w in (417) is calculated from 



in which Ax is the ^integration step size in the current step and Axq that in the 
preceding step. Thus, w never exceeds 0.5, and equation (417) provides some 
smoothing of fluctuations in (dS*/dx) c . if the current step size Ax is much 
smaller than the value Axq for the preceding integration step, which can happen 
if the current flow point is a model point or if the numerical integration of the 
flow equations is encountering severe convergence problems, then the weight w 
applied to the current (dS*/dx) c value is very small. This algorithm avoids 
difficulties which could result from loss of accuracy in the subtractions indi- 
cated in (416) wheh x and x Q are very nearly equal. 

The above-described procedure for switching on the coupling between the 
boundary layer and the inviscid flow is found to work reliably in flow calculations 
based on conventional chemical non-equilibrium models. However, when attempts are 
made to compute flow solutions including the boundary layer for gases represented 
by the electronic non-equilibrium model, the solutions usually fail at the point 
where the coupling is switched on. Such failures occur because the disturbance 
produced by increasing dS*/dx in accordance with (417) and (418), small though 
it is, upsets the precarious stability of the electronic non-equilibrium model. 

Boundary layer calculations are performed both at the intermediate points in 
the Treanor-Runge-Kutta integration step and at the end of each successful step. 
Now, the method used in the boundary layer calculations involves the evaluation 
of a definite integral, equation (170) or (172), in which the variable of inte- 
gration is the streamwise coordinate f lying along the nozzle surface. Also, f 
itself is evaluated as a definite integral. To avoid contaminating these inte- 
grals with inaccurate data from unsuccessful integration steps and intermediate 
points in successful steps, the "permanent" values of the integrals are incre- 
mented only following successful full steps. At intermediate points in a step, 
the integrals are incremented on a temporary basis for use in computing the 
boundary layer displacement thickness, but are restored to their initial values. 
Also, d8* / d x is not computed at intermediate points. 


The calculation of the displacement thickness 5* (Sections 5.5 and 5.6) 
involves the quantity dKrvM/dx, in which M denotes the Mach number. During the 
non-equilibrium solution, this derivative is evaluated as follows: 

d In M 1 du 1 da 1 du 1 dT (420) 

dx u dx a dx u dx 2T dx 


where the second form is based on the approximation a « t for the sound 
speed a. The velocity derivative is obtained from equation (245): 


du 

dx 



(421) 
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For chemical non-equilibrium gas models, dh Q /dx = 0 and dh/dx is given by (291). 
For electronic non-equilibrium models, dhg/dx is given by (333) and, from (334), 


dh 

aT = } ' e c p e 



(422) 


Thus, dEn M/ dx is evaluated from the current values of dT/dx, dyj/dx , and, in the 
case of electronic non-equilibrium models, dT e /dx and dh 0 /dx. These derivatives, 
in turn, are obtained by solving equations (293) and (299) (or (328), (331), (332), 
and (333) for electronic non-equilibrium models). These equations for the deriva- 
tives involve A e and d ( n kj dx , which depend on S* and dS*7dx. Thus, the calcula- 
tion of S* at any flow point requires knowledge of the value of 8* at the same 
flow point. This impasse is broken by using an iterative self-consistent solution 
for 8* and the derivatives dT/dx, dyj/dx , dT e /dx, dhQ/dx D f the flow variables. 

The convergence criterion in the iteration is that the newly computed 8* differ by 
no more than 1 percent from the previous value. Whenever calculations of the 
derivatives are required, in a non-equilibrium solution including the boundary 
layer, the derivatives are first computed assuming the 8* value 8 p * left in storage 
by the last previous calculation. Then a new value 8j* for the thickness is com- 
puted by calling the boundary layer routine, and is compared with the previous 
one. If Si* and S p * differ by no more than 1 percent, the new value is accepted 
and used. If they differ by more than this amount, the entire calculation is re- 
peated, using Si in the rate calculation, to obtain a further improved estimate 
S 2 *. If the convergence test is satisfied, S 2 * is accepted as the displacement 
thickness. If not, the following rapidly converging algorithm is used to obtain 
a more accurate estimate of S*. Denote the input value of displacement thickness 
assumed in the rate calculation by 8 -* , and the output value obtained from the 
boundary layer routine by S Q *. Then the calculational procedure defines a func- 
tional relationship 8 0 * = f(Sj*) . The condition for a self-consistent solution is 
that 8 0 * = 8f . The functional relation is approximated by a straight line passing 
through the two points already computed (or the two most recent points if the 
algorithm is applied more than once) . The equation of the straight line is 


V - S o2* 


S o2 - Sol 


Sf - 8; 


i2 


S]2 - S; !* 


(423) 


The intersection of (423) with 
at 8 0 * = Sj* = 8* , where 

* s c ; Sji - s 0l * s i2 * 


«il - S i2 


- Sol +s o2 


the self-consistent solution line is found to be 


(424) 


Figure 27 illustrates the geometry of this solution. The first application of 
(424) almost always gives a S* value which satisfies the self-consistency con- 
dition to within 1 percent. The program allows a maximum of four iterations 
using (424) . 
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7.7 Example - Nf^ Recombination 


To illustrate the non-equilibrium flow calculations performed by NATA, the 
code has been used to compute the flow and concentration changes in two of Wegener's 
classic experiments (refs. 46, 47) on the three-body recombination of nitrogen di- 
oxide , 

2 N0 2 + N 2 F=*N 2 0 4 + N 2 (425) 

Wegener's experiments were performed in a small wind tunnel of rectangular cross 
section. The flow geometry was modeled in NATA as a two-dimensional nozzle with 
the profile shown in figure 28. The flow in this tunnel was found to be one- 
dimensional to within the accuracy of the measurements. The experiments analyzed 
were those designated C and F in reference 46. The conditions in these experi- 
ments were as summarized below: 


Experiment 

< V >0 

Po 

T o 


2 

percent 

atm 

°K 

C 

0.0100 

2.00 

400 

F 

0.0500 

2.16 

402 


Of the three gas species appearing in (425) only N2 is a standard species in 
NATA. Thermo fits for NO2 and N2O4 were obtained by fitting the data in the JANAF 
tables (ref. 48) over the temperature range 200° to 400° K. The coefficients ob- 
tained for use in equations (33) and (34) were as follows: 


Coefficient 

N0 2 

N 2 °4 

a 

4.003 

3.553 

b 

-3.75 x 10~ 4 

1.1625 x 10 -2 

c 

2.45 x 10 -6 

-4.55 x 10 -6 

d 

0 

0 

e 

0 

0 

k 

5.945 

10.028 


Hq° (kcal/mole) 8.586 4.473 

The rate constant for the reaction (425) was taken to be 3 x 10^ 4 cm*Vmole 2 -sec 
(ref. 47). 

' Frozen, equilibrium, and non-equilibrium flow calculations were performed 
for each of the two experiments. The results are compared with Wegener's data 
in figures 29 and 30. The abscissa in these figures is distance along the 
nozzle. The ordinate is the ratio [N02]/[N02] 0 of the N °2 concentration to its 
value in the upstream reservoir. This was calculated from quantities available 
in the NATA output as 
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Figure 29 COMPARISON OF NATA RESULTS WITH WEGENER'S DATA ON N0 2 RECOMBINATION 

(EXPERIMENT C) 
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Figure 30 COMPARISON OF NATA RESULTS WITH WEGENER'S DATA ON NO2 RECOMBINATION 

{EXPERIMENT F) 




p 


0 


( 426 ) 


[no 2 ] 

[no 2 ] 


x no 2 w 
. 0 


The curves in the figures give the NATA results. The points are Wegener's data 
(ref. 46), which are based on measurements of optical transmission at wavelengths 
of 3950 to 4750 A. In this wavelength range, N 2 and N 2 O 4 . are transparent, while 
NO 2 absorbs. The figures show that the NATA non-equilibrium flow calculations 
are in close agreement with the experimental data in both cases. In experiment 
F, the flow is nearly in equilibrium throughout the tunnel, but a major departure 
from equilibrium occurs in experiment C. 
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8. CONDITIONS ON MODELS 


The NATA code provides calculations of the heat flux and pressure on models 
inserted into the free-stream flow at specified locations. The code treats two 
general types of model configurations : the stagnation point of an axisymmetric 

or two-dimensional probe, and the flat surface of a blunt wedge. 

The points in the free-stream flow at which model condition calculations 
are to be done may be specified in either or both of two ways: 

(1) A list of up to 20 test section diameters (TSDIAM) may be read in. 

Model condition calculations are done at the positions (downstream of 
the throat) at which the nozzle diameter has the values thus specified. 

(2) An initial axial coordinate value (XM0DP1) , a final value (CXMAXI), and 
a number of model points (NM0DPT) may be read in. Model condition 
calculations are then done at the specified number of points, including 
the initial and final positions and a sequence of intermediate points 
whose axial coordinates form a geometric progression. For NM0DPT =1, 
calculations are done only at the initial point. 

The influence of the test models upon the free-stream flow is not considered. 
Thus, when calculations are done for several model locations, each model is as- 
sinned to be inserted into the previously undisturbed free stream. 


8. 1 Stagnation-Point Calculations 

The calculation of stagnation-point conditions begins with an approximate 
normal shock solution. The heat flux and shock standoff distance are then cal- 
culated for both flat-faced and hemispherical axisymmetric models. Optionally, 
the calculations can be done instead for two-dimensional models with flat and 
cylindrical leading edges. These calculations are carried out assuming either 
equilibrium flow behind the shock or a frozen shock with the species mole frac- 
tions behind the shock equal to those ahead of it. Optionally, both types of 
normal shock calculation can be done. 

Section 8.1.1 discusses the normal shock solutions. Section 8.1.2 the heat 
flux calculations, and Section 8.1.3 the stagnation-point velocity gradient. 
Section 8.1.4 specifies the model heat flux outputs provided, and Section 8.1.5 
discusses low-density effects which limit the validity and accuracy of the 
stagnation-point heat flux calculations. 


8.1.1 Normal Shock Solutions 

The main purpose of the stagnation-point model-condition calculations is to 
predict the stagnation point heat flux and the stagnation pressure on the model. 
The objective of the normal shock solution is to determine the conditions at the 
stagnation point of the inviscid flow, just outside the boundary layer on the 
model. The free stream conditions ahead ef the shock will be designated by 
subscript 1 , conditions just behind the shock by subscript 2 , and stagnation 
conditions by subscript s . 
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The conservation relations for the normal shock are: 
Pi + Pi“i^ — P 2 + P2 u 2^ 

Pi U 1 = P2 u 2 

1 2 1 2 
hi + jui = h 2 + y u 2 


(427a) 

(427b) 

(427c) 


The conditions ahead of the shock are known from the nozzle flow solution. The 
conditions at the inviscid stagnation point are assumed to be related to those 
just behind the shock by 

h s = h 2 + J u 2 2 (428a) 


Ps 


P2 



(428b) 


The first of these relations (428) is exact. The incompressible Bernoulli equa- 
tion (428b) is only an approximation to the actual' compressible isentropic flow 
relation, but is reasonably accurate as long as the shock Mach number is not too 
low. In NATA, model condition calculations are not done for points at which the 
free-stream Mach number is less than 1.5. At this Mach number, for a perfect 
gas with y = 1.4 (for example), equation (428b) is accurate to within about 3.5 
percent. For higher free-stream Mach numbers, the accuracy is still better be- 
cause the Mach number behind the shock is lower. For example, at Mach 5 the 
error is only 0.5 percent for a perfect gas with y = 1.4. 


The gas equation of state has the form 

P = p(T.p) (429a) 

h = h (T, p) (429b) 


In the case of a frozen shock, the relations (429) are assumed to be the equa- 
tions for a thermally perfect but calorieally imperfect gas: 



(430a) 


h = h(T) (430b) 

For the equilibrium shock, the relations (429) are embodied in a subroutine 
(EQCALC) which computes the thermochemical equilibrium state of a gas mixture at 
specified temperature and pressure. In this case, neither of the equation of 
state relations can be expressed in analytical form. 

NATA calculations of the stagnation conditions on a model are carried out in 
two stages. First, the conditions immediately behind the shock are computed by 
simultaneous solution of the Rankine-Hugoniot conservation relations (427) and 
the equation of state (429) . Then the -stagnation conditions are calculated by 
simultaneous solution of equations (428) and (429). 
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The equations (427) for the normal shock can be rewritten in terms of the 
density ratio 

e P ~ P1 ; P2 (431) 

in the following form: 

P2 = Pi + Pi u l 2 ( 1 ~ e p^ (432a) 

h 2 = h x + i Ul 2 (l- v 2 ) (432b) 

In (432) , equation (427b) has been used to eliminate the velocity u 2 from the 
remaining two equations (427a), (427c). 

For both the equilibrium shock and the frozen shock, the solution involves 
an iteration on the temperature T 2 of the gas immediately behind the normal shock. 
For each estimate of T 2 , the static enthalpy h2 is calculated both from (432b) 
and from the equation of state relation (429b) or (430b) . The correct temperature 
T 2 is the one for which these two enthalpy values are equal. 


In the case of a frozen shock, the gas enthalpy (430b) depends on tempera- 
ture only, not on pressure. It is computed from the species enthalpies Hj(T 2 ) 
using equation (275). In the solution for the frozen shock, the mole fractions 
are assumed constant across the shock, so that the free stream values are used. 
To calculate the second h 2 value from (432b) , the density ratio fp must be known. 
In the case of the frozen shock, since the gas is assumed thermally perfect and 
since the molecular weight W behind the shock is the same as that ahead of the 
shock, the density P 2 behind the shock is given by 


P2 = Pi 


P2 T l 

Pi t 2 


(433) 


From (431) and (433) 

Pi t 2 

V = ^ ' ~1 


(434) 


Substitution of (434) into (432a) gives a quadratic equation for the pressure P2 : 


p 2 2 _ p T p 2 + Ap = 0 

where 


p T = Pi + Pi Uj 2 

Pi “i 2 Pi t 2 



The applicable solution of (435) is 
1 

P2 = 7 


p T 




Pq- 


-4 A. 


(435) 

(436a) 

(436b) 


(437) 
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For each estimate of the temperature T 2 in the iterative solution of the equa- 
tions for the frozen normal shock, P 2 is calculated from (437) and e p is then 
obtained from (434) . This *p value is then used to compute h2 from (432b) . The 
resulting h 2 value is compared with the one from (275) to determine whether con- 
vergence has been achieved and if not, to select an improved estimate of T 2 for 
the next iteration. 

In the case of an equilibrium shock, the solution is slightly more compli- 
cated because the enthalpy (429b) depends upon the pressure P 2 , and because the 
calculation of ( p for each estimate of T 2 cannot be performed analytically but 
requires an inner iteration. For each trial value of T 2 , the corresponding 
pressure P 2 and density ratio e p are obtained by solving equations (429) , (431) 
and (432a) , as follows : 

1. An initial estimate of ( p is chosen and is used to calculate an initial 
estimate of P2 from (432a) . 

2. The equilibrium equation of state (429a) is used to compute 
p2 = p (T 2 , p 2 ) 

3. The new density ratio 

e p = PVP2 

is calculated. 

4. The new ratio e p is compared with the old value e p to determine whether 
the inner iteration has converged. The criterion is that the two values 
be equal to within 0.1 percent. 

5. If not, the quantity 
a = p 2 /P2 


is computed, and a new estimate of p 2 is obtained from 


p 2 = 




(438) 


Equation (438) is equivalent to (437) with 

Ap = p 1 2 u 1 2 a 

Thus, (438) is an approximation to the pressure p 2 corresponding to the 
assumed temperature T 2 , based on neglecting the pressure dependence of 
the gas composition. 

6. The steps beginning with (2) are repeated using the new P 2 from step 
(5). 
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* 

This inner iteration to obtain ep and P2 is repeated for each step of the main 
iteration on T 2 . The logic of the iteration to determine T 2 is the same as for 
the frozen shock. 

Once the conditions P2 , T 2 , p 2 behind the shock have been obtained, the 
conditions at the inviscid stagnation point on the model are calculated by solv- 
ing equations (428) and (429) . The stagnation pressure is calculated directly 
from (428b) , which is rewritten in the form 
1 2 

p s = p 2 + j Pl u 1 e p (439) 

The stagnation temperature T s is then obtained by an iterative solution of (428a) 
and (429b) or (430b) . 


8.1.2 Stagnation-Point Heat Flux 


NATA calculations of stagnation-point heat transfer to models are based on 
a correlation formula which incorporates results of analyses by Fay and Riddell 
(ref. 49), Bade (ref. 50), Finson and Kemp (ref. 51) and Goulard (whose work is 

52, pp. 92-93). Fay and Riddell's correlation 


reviewed by Dorrance, ref. 
%u,w 




0.763 N Pr l 


0.4 


e, w 


Pe 

PwPw 


0.4 


1 + (N 


f 


Le 


1 ) 


n D 


(440) 


(ref. 49) is widely used for calculating stagnation point heat transfer in high- 
temperature air. The non-dimensional groups appearing in this formula are de- 
fined as follows: 


Nusselt number based on property values at the wall 
<ls x c pw 


N 


Nu, w ~ 


K w (h e - h w ) 


(441) 


Reynolds number based on property values at the wall 


N 


Re, w 


Pw u e x 
Pw 


(442) 


Frozen Prandtl number (assumed constant) 
c p P 

N Pr = v ' 


Lewis 


N 


Le = 


number (assumed constant) 
D 12 P c p 
K 


(443) 


(444) 
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The other notations used in (440) to (444) are: 

Cp = specific heat (chemically frozen) 

D 12 = binary diffusion coefficient for atoms and molecules 

f = 0.52 for an equilibrium boundary layer, 0.63 for a frozen one 

h = enthalpy 

h D = enthalpy of dissociation per unit mass for the gas in the external 
flow 

K = thermal conductivity (chemically frozen) 
q s = stagnation point heat flux 

u = velocity component in the direction parallel to the body surface 

x = radial distance along the body surface from the stagnation point 

ft = viscosity 

P = density 

The subscript notation is 

e = conditions at the stagnation point of the external flow 


w = conditions at the body surface 


It is easy to show from the above definitions that 



(445) 


Thus, for constant Prandtl number the Fay-Riddell formula (440) can be rewritten 
in terms of Nusselt and Reynolds numbers based on property values at the edge of 


the boundary layer in the form 

N Nu, e o A / Pw ^ 

-== 0.763 N Pr °-4 


l + (N Le f -1) 


The Fay-Riddell' correlation (440) or (446) is a curvefit to the results of 
49 numerical solutions of the stagnation-point laminar boundary layer equations. 
The solutions included cases in which the boundary layer was assumed to be in 
chemical equilibrium, cases in which finite rates were used for the dissociation/ 
recombination reactions, and cases with a chemically frozen boundary layer. In 
all cases, the viscosity was assumed to be given by a Sutherland formula. 




( 447 ) 
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with po = 1.72 x 10“4 poise, Tq = 273° K, and Tg = 113° K, and the Prandtl num- 
ber was assumed to have the constant value 0.71. The Fay-Riddell formula (440) 
or (446), with these gas transport properties, with f = 0.52 (equilibrium 
boundary layer), and with a Lewis number NL e = 1.4, has been found to correlate 
measurements of stagnation-point heat transfer in air to within the experimental 
accuracy of about ±20 percent (ref. 52, p. 91; fef. 53), at temperatures up to 
the onset of ionization effects. 

Although the Fay-Riddell correlation was developed specifically for disso- 
ciated air, it can also be used to estimate stagnation-point heat transfer in 
other gas mixtures. This is done in NATA. The transport properties p , K , Np r , 
^Le» appearing in (446) and in the defining equations (441) and (442) for the 
Nusselt and Reynolds numbers, are computed by the NATA transport subroutine 
package (Section 3) . 

Air is treated in the same way as other gases; i.e., transport properties 
computed by 'the code are used in place of the properties, such as (447), assumed 
by Fay and Riddell in their calculations. For this reason, a minor modification 
of the correlation formula (446) is required. The Sutherland formula (447) for 
the viscosity of air is quite inaccurate at high temperatures. For example, at 
7000° K and 1 atm pressure, the value based on (447) is 35 percent lower than 
the value given by Yos' calculations (ref. 54). At high temperatures, the 
Sutherland formula is practically equivalent to a power-law viscosity formula, 
p = AT i, ) with the exponent <u = 0.5 , whereas realistic calculations for equilib- 
rium air (ref. 54) can be approximated by the power law with co ~ 0.7 in the range 
from room temperature up to 5000° K. At higher temperatures, the viscosity first 
increases still more rapidly, then passes through a maximum and falls off. 


The value of the exponent <o affects the dependence of the stagnation-point 
heat flux upon the enthalpy ratio g w = h w /h e across the boundary layer (ref. 50). 
For co = 0.5, the correlation parameter NNu,e/V N R e , e decreases with increasing 
enthalpy ratio. This is the dependence which Fay and Riddell fitted using the 
factor (p w n-v/pe p e ) 0-1 in (446). However, their representation of this dependence 
is valid only for <u = 0.5 (or for the Sutherland law). For co = 1, (p w p^/ p,, p„) 
is equal to unity for all values of the enthalpy ratio, but 1% U) e /\/NR e " e increases 
slightly with increasing g w . For co - 0.7 to 0.8, (p w p,,/^ p e ) varies with g w but 
Nnu, e/V N Re, e is practically independent of g w . 


Since the viscosity of air, as calculated in NATA, is approximated by a power 
law with co — 0.7 for temperatures up to 5000° K, the results of reference 50 indi- 
cate that the dependence of the correlation parameter N Nu g/V^Re e u P on the en- 
thalpy ratio should be weak. For this reason, the Fay-Riddell correlation (446) 
is adjusted for the use of the real air viscosity in place of the Sutherland law 
(447) by omitting the factor (p w p,/Pe representing the dependence on the 

enthalpy ratio across the boundary layer: 


N, 


Nu, e 


Ni 


Nu 


\f~N 


Re, e 




1 + (N Le f _ 1) 


n D 


' o L 


where 


/ n Nu \ 

\ V N R e / o 


c <i) N Pr,w' 


0.4 


(448) 


(449) 
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is the laminar heat transfer parameter for a low-temperature flow with constant 
Prandtl number and pp - constant through the boundary layer. In (449), C(j) is 
a coefficient which depends upon body geometry. For a stagnation point in two- 
dimensional flow (j=0 ) , according to Squire (ref. 55), 

C(0) = 0.570 (450a) 
For an axisymmetric stagnation point (j=l), from Sibulkin (ref. 56), 


C(l) = 0.763 (450b) 

These values (450) can also be obtained from the results of Dewey and Gross 
(ref. 26).* 


Equation (448) shows that the laminar heat transfer parameter, evaluated 
using gas properties in the external flow, may be approximated as the product of 
the low- temperature correlation function (449) and a correction factor repre- 
senting the effects of heat transfer by diffusion and recombination of atoms. 
Substitution of the definitions for Nfcj u , e and Nr 6i e (analogous to equations 
(441) and (442) for N]sj u , w and NR e , w ) into (448) gives the following expression 
for the stagnation-point heat flux: 


% = 


N 


C(j) 

(L6 


Pr 



/ duA 

Pe Pc 

/ e \ 

\ dx ) s 


~il/2 


(h e -h w ) 


1 + <N L , 


f 


1) 


(451) 


In the vicinity of the stagnation point, the velocity u e is proportional to the 
radial distance x from the stagnation point. Thus, in (451) the ratio Ug/x has 
been replaced by the stagnation-point velocity gradient, ( dug/dxlg 


A further minor modification of the Fay-Riddell formula has been made to 
permit inclusion of the effects of surface catalytic efficiency in the calcula- 
tions. In the case of a frozen boundary layer, (446) and (451) assume that the 
model surface is perfectly catalytic for recombination of atoms. In NATA, (451) 
has been modified to allow for surface catalytic efficiency less than 1, using 
an approach suggested by Goulard (ref. 52, pp. 92-93): 


= 


C(j) 
0j$ 


N 


Pr 




1 + (<S N Le f - 



<h e - h w> 


(452) 


Here 4 is a catalytic factor which can take values between 0 (no catalysis) and 
an upper limit which is less than 1 for a perfectly catalytic surface. Equation 
(452) is used only for the frozen boundary layer. For the equilibrium boundary 
layer, the factor $ is omitted because, in that case, recombination of atoms is 
accomplished by reactions in the boundary layer regardless of the properties of 
the model surface. For $ = 0, (452) makes the heat flux proportional to 
(h e — hj} — h w ) . According to Goulard's analysis, $ is given by 

(N L e/N Pl ) 2/3 [2(du e /dx) sPeMe ] 1/2 


Pw ^Rw 




*T o get (450) fro m the data in De wey and Gross, it is necessary to note that the Falkner-Skan parameter p is 1 for j =0 and 0.5 for j = 1 , 
and that N Nu /\/N Re = vTTi 6' (0) for 1^ = 1. 
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where kp w is the surface recombination rate constant, defined such that the mass 
of atoms recombining on the surface per unit area per unit time is kRw times the 
concentration of atoms in mass per unit volume in the gas at the surface. The 
units of k Rw are cm/sec. 

In NATA, the dissociation enthalpy is calculated as 


h ° = T ^2 Xi Hi0 ° (454) 

atoms 

in which w is the mean molecular weight, Xj is the mole fraction for the ith 
species, and H;o° is the molar enthalpy of formation for the species. Only neu- 
tral atoms are- included in the sum. The Lewis number Nl c is calculated for an 
atom-molecule diffusion process which is assumed to be dominant or characteristic; 
in the standard air models, N P e is computed based on the O-N 2 diffusion coeffi- 
cient. Thus, the significant recombination process is assumed to be 20— ►O 2 . 

In the planetary atmosphere (C02-N2“Ar) models, the bracketed factor containing 
the Lewis number in (448) , (451) , and (452) is omitted since there is no theo- 
retical support for applying it to a gas mixture so different from air. 

The boundary layer solutions of Fay and Riddell did not include effects of 
ionization of the air. Flows with copious free-stream ionization are within the 
intended scope of application of the NATA code. Gas ionization affects the 
structure and properties of a laminar boundary layer primarily through its effects 
on the gas transport properties. The main changes in transport properties due to 
ionization are as follows: 

(1) The chemically frozen thermal conductivity Kf increases greatly as a 
result of the high conductivity of the free electrons. 

(2) Energy transport by ambipolar diffusion and recombination of electron- 
ion pairs occurs. This process further increases the effective thermal 
conductivity. However, this reaction conductivity is considerably less 
important than the corresponding effect due to dissociation and recom- 
bination of diatomic molecules in air, because the large charge exchange 
cross sections impede the diffusion. 

(3) The viscosity decreases as a result of the large momentum-transfer 
cross section for collisions of charged particles. 

As a result of these changes, the (frozen) Prandtl number Np r = c p f p/Kf de- 
creases markedly as the amount of ionization increases. Thus, when the external 
flow is ionized to a substantial extent, the frozen Prandtl number is considerably 
smaller at the outer edge of the boundary layer than at the cold surface, where 
the ionization is negligible. This variation of the frozen Prandtl number is a 
new feature, not previously encountered in studies of lower temperature boundary 
layers. In air, for example, the equilibrium Prandtl number ( c p) e q/ i / K eq (where 
K e q includes the reaction conductivity) has peaks and dips due to the reactions, 
but the frozen Prandtl number varies by only a few percent in the temperature 
region below the onset of ionization. 
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Available stagnation point boundary layer solutions have been studied to 
obtain a correlation of the effects of free-stream ionization upon stagnation- 
point heat transfer. These solutions include the calculations of Finson and 
Kemp (refs. 51, 57) for the noble gases, those of DeRienzo and Pallone (ref. 58) 
for air, and those of Fay and Kemp (ref. 59) for nitrogen. 

Finson and Kemp (refs. 51, 57) computed solutions of the laminar stagnation- 
point boundary layer equations in axisymmetric and two-dimensional flows of 
ionized argon and xenon. They included energy transport by ambipolar diffusion 
and recombination, and obtained two sets of solutions, one based on local thermo- 
chemical equilibrium through the boundary layer, the other based on frozen 
chemistry but assuming a perfectly catalytic wall. The cases computed were 
chosen to simulate situations for which experimental data were available in the 
literature. 

Figure 31 presents a correlation of Finson and Kemp’s results for the equi- 
librium boundary layer. The abscissa in this figure is the ratio of the frozen 
Prandtl number across the boundary layer. In all cases, the wall temperature 
was 300° K; thus, N Pfi w was constant at its low- temperature value of 2/3 for 
monatomic gases. The frozen Prandtl number Np r> e in the external flow was 
calculated from the transport property formulas assumed in Finson and Kemp's cal- 
culations, rather than from other, possibly more accurate data, because the ob- 
ject of the study was to correlate the calculated heat flux with the transport 
properties upon which it was based. 


The ordinate in figu re 3 1 is the ratio of the non-dimensional laminar heat 
transfer parameter, ( N Nu /\/N Re ) e , based on gas properties in the external flow, 

to the value ( n Nu/ V N R e ^0 for a low- temperature flow, given by (449). Bade (ref. 


50) has shown that this ratio is approximately equal to unity for stagnation 
point boundary layers in helium, neon, argon, krypton, and xenon at temperatures 
up to the onset of free-stream ionization. Figure 31 shows that, in the ioniza- 
tion region, this ratio varies over a wide range. The straight line drawn to 
fit the data points in the figure represents the equation 


(%/ VNjTe)e / N Pr, e 

(n Nu/V^rPo \ Npr - w 


(455) 


Substitution of (441), (442), and (449) into (455) gives the following formula 
for stagnation point heat transfer in a noble gas, valid from low temperatures 
up to and including the ionization region: 


'J P e 4 e (d u e/dx) s 

% = C(j) : (h e -h w ) (456) 

(Npr,wNp r , e >°- 3 


This formula is used in NATA for calculations of stagnation point heat transfer 
in helium and argon. Since the effect of ionization is to decrease N Pl) e , 
equation (456) shows that free-stream ionization increases the stagnation point 
heat flux above the value predicted by the theory of reference 50 for a gas of 
neutral atoms at the same enthalpy.. 
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= (N Nu /V^ e )e/(NNu/VNRe)0 



84-1802 N Pr,e/NPr,w 


Figure 31 CORRELATION OF STAGNATION-POINT HEAT TRANSFER IN 

IONIZED NOBLE GASES 




Theoretical calculations of stagnation point heat transfer are also avail- 
able for partially ionized air and nitrogen. DeRienzo and Pallone (ref. 58) 
computed the heat flux through an equilibrium boundary layer in air. They took 
energy transport by diffusion and reactions into account by using the equilib- 
rium Prandtl number and omitting the diffusion equations from the system to be 
solved. Their results for equilibrium air with no blowing at the wall are 
represented by the circles in figure 32. Fay and Kemp (ref. 59) performed cal- 
culations for nitrogen, using a model in which the reaction conductivity due to 
the ionization reaction Nf^N^+ewas neglected. This approximation was justified 
by reference to the large charge-exchange cross section for nitrogen. Fay and 
Kemp explicitly included energy transport by diffusion of neutral atoms and their 
recombination to form diatomic molecules. Their results are represented by the 
diamonds and triangles in figure 32. 

Both reference 58 and reference 59 tabulate their results in terms of the 
parameter (NNu/V N R e )w » based on gas properties at the cold wall. For the 
preparation of figure 32, these data were "converted into ( N N u /V N R e )e values with 
the aid of equation (445) . In each case, the transport properties assumed in 
the reported boundary layer solutions were used in equation (445) and in evaluat- 
ing the abscissa of figure 32, since the object of the study was to correlate 
calculated heat transfer data with the transport properties upon which the calcu- 
lations were based. In the case of reference 58, the surviving documentation is 
not sufficient to specify the assumed transport properties exactly, so that there 
is some uncertainty as to where the points should lie. Fay and Kemp provided 
enough information to permit reconstructing the transport property values for 
each of their solutions. 

The scatter of the points in figure 32 is considerably greater than that in 
figure 31, but the overall trend is similar. The straight line in figure 32 
again represents equation (455) . The fit would be a little better if the expo- 
nent of 0.7 in (455) were replaced by a slightly lower value, e.g., 0.65. How- 
ever, (455) accounts for the main part of the variation of the ordinate. The 
uncorrelated variations probably result mainly from energy transport by diffusion 
and recombination of atoms. Trend lines drawn through the Fay-Kemp results for 
L = 0.3 and L = 1 (triangles and inverted triangles) bracket nearly all of the 
data from both references 58 and 59. The parameter L is not the Lewis number as 
defined by (444) , but is a related quantity which (unlike Nl c ) is approximately 
constant over the temperature range in which nitrogen molecules- dissociate. Fay 
and Kemp's results indicate that the stagnation point heat flux varies roughly 
as L 


In NATA, stagnation point heat transfer in air and other molecular gases 
under conditions of free stream ionization is calculated by combining the 
ionization effect represented by equation (455) with the atom diffusion and re- 
combination effect given by (448) (and adjusted for surface catalytic efficiency 
in (452) ) . This combination of correlations gives 


(N N u/V%e>e _ / N Pr, e 

( N NuA/NRe>0 \ Npr > w 


1 + (4> N Le f 



(457) 
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= (NNuA/NR^) e /(N|\j u /\/NRe)o 



84-1803 


Npr,e^Pr,w 


Figure 32 CORRELATION OF STAGNATION-POINT HEAT TRANSFER IN 
IONIZED AIR AND NITROGEN 
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The corresponding stagnation point heat flux formula is 


c (j) l du e\ 

(N N >0.3 J Pet * e l dx J 
(N Pr, w N Pr, e> yl \ ' ' 


1 + (* N Le f - 1) 


a D 


h e h w 


(h e -h w ) ( 458 ) 


This differs from equation (456), which is used for the noble gases, only by the 
factor in square brackets representing the effects of atom diffusion and re- 
combination. 

'"J- 


8.1.3 Stagnation Point Velocity Gradient 

The stagnation-point velocity gradient (du e /dx) s appearing in (452), (456), 
and (458) depends upon the free-stream velocity uj , the model nose radius % , 
the shock density ratio e p , and the model shape. It also depends weakly upon the 
gas equation of state, i.e. , the effective specific heat ratio, y . 


For hemispherical models in axisymmetric flow and cylindrical models in two- 
dimensional flow, NATA computes the velocity 'gradient using the modified 
Newtonian approximation (ref. 60) : 


P-Pl , 

= cos"^ /3 

Ps - Pi 


(459) 


where p is the local pressure on the body surface at a point where the incident 
flow ahead of the shock makes an angle /3 with the normal to the surface. As 
before, pi is the free-stream static pressure ahead of the shock and p s the stag- 
nation pressure on the body. Near the stagnation point, the pressure and flow 
velocity are related by the incompressible Bernoulli equation 


P s = P + 



(460) 


and the distance along the surface from the stagnation point is given by 


x = R n /3 - R„ sin /3 

where R n is the body radius of curvature at the stagnation point, 
of (459), (460) and (461) gives 


V dx )s R n i Ps V Ps / 


(461) 
Combination 

(462) 


Equation (459) is an empirical pressure distribution which is reasonably 
accurate for spheres. Figure 33 compares the non-dimensional velocity gradient 
parameter. 


Rn 

U 1 


du^ 


dx 



(463) 
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based on the modified Newtonian formula (462), with experimental data from Boison 
and Curtiss (ref. 60), Crawford and McCauley (ref. 61), and Reinecke (ref. 62).* 
The figure also shows results of unpublished calculations carried out with the 
Mbretti-Bleich blunt-body computer program (ref. 63) by R. C. Boger of Avco 
Systems Division. The curve represents equation (463) . The modified Newtonian 
formula agrees with the experimental data to within their scatter, and lies less 
than eight percent below the trend of the theoretical results from the blunt-body 
inviscid flow program. Since the stagnation point heat flux depends upon the 
velocity gradient as (du e /dx) s 1/2 , the accuracy of (463) for models with spherical 
noses appears adequate. 

Information surveyed by Hayes and Probstein (ref. 64, pp. 422-423) suggests 
that the pressure distribution in two-dimensional hypersonic flow over a circular 
cylinder is given approximately by equation (459) . The velocity gradient param- 
eter ( R n / U i) (du e /dx) s is apparently a little lower for cylinders than for spheres. 
This difference, if real, is neglected in NATA. 

Figure 34 summarizes some available data on the stagnation point velocity 
gradient in axisymmetric flow over a flat-faced cylinder. The experimental data 
from reference 60 lie about 20 percent below the results of Vinokur's (refs. 65, 
66) approximate constant density solution. The unbroken curve shown represents 
the analytical relation 

5: . S (464> 

U 1 \ dx / s 1.34+ 4.2 (v^-0.4) 2 

which is a curvefit to Vinokur's results for the flat-faced cylinder, divided by 
a factor of 1.2 to improve the fit to Boison and Curtiss's experimental data. 

In (464), Rb denotes the radius of the circular flat face. 


Two-dimensional flow over a flat-faced plate has been analyzed by Cole and 
Brainerd (ref. 67) using the constant density approximation. According to Hayes 
and Probstein (ref. 64, p. 310), Cole and Brainerd's calculations give the 
stagnation point velocity gradient in the form 



where R n is the half width of the flat face. Hayes and Probstein (ref. 64, pp. 
426-428) have compared Cole and Brainerd's results with those from other 
calculations . 


Since the body radius Rb is equal to the nose radius of curvature R n for 
spherical and cylindrical models, equations (463), (464), and (465) all indicate 
that the stagnation point velocity gradient (du e /dx) s is equal to ( U l/Rt>) times 
a non-dimensional function of the inviscid flow parameters on the stagnation 
streamline. Thus, the stagnation point heat flux (458) can be expressed in the 
form 


(b e h w ) (466) 


= 


C(j) 

(N p Np g) 1 


0.3 


'Jp e l l e' 1 l G 


1 + ($ N Le f - 1) 


“D 


*The authors are indebted to W. G. Reinecke of Avco Systems Division for assembling the data plotted in the figure. 
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EXPERIMENT (BOISON AND CURTISS) (ref. 60) 
VINOKUR (CONSTANT-DENSITY) (refs. 65, 66) 



FLOW OF AIR OVER A FLAT-FACEL^YLINDER 




( 467 ) 


where G denotes the non-dimensional velocity gradient parameter: 



If the catalytic factor $ is considered to be a given constant, then the right- 
hand side of (466) is independent of body radius. However, if $ were assumed to 
be given by equation (453) , the right-hand side of (466) would depend explicitly 
on R-b . The present version of NATA treats $ as an input constant, so that (466) 
is independent of body radius. 

8.1.4 5 Model Heat Flux Outputs 

.J 

The stagnation point heat flux on the model is printed out as q s , for 
up to eight cases. The following four values are printed when equilibrium 
shock calculations are done: 

[q s VR- b ] EFE Equilibrium shock, flat-faced model, equilibrium boundary 
layer 

[q s V^b^EFF Equilibrium shock, flat-faced model, frozen boundary layer 

KVR^EHE Equilibrium shock, hemispherical model, equilibrium 
boundary layer 

[q s V^b'JEHF Equilibrium shock, hemispherical model, frozen boundary 
layer 

If frozen shock calculations are requested, the following four values are printed: 

[q s V RJfFE Frozen shock, flat-faced model, equilibrium boundary layer 

[q s V R b^FFF Frozen shock, flat-faced model, frozen boundary layer 

[q s V^blpHE Frozen shock, hemispherical model, equilibrium boundary 
layer 

[q s V R b^FHF Frozen shock, hemispherical model, frozen boundary layer 

In addition to these values, NATA provides two other heat fluxes at each 
model point. One is based on a simple formula proposed by Hiester of Stanford 
Research Institute (ref. 68): 

W0M = °- 0417 ^ ^ " h w> (468) 

where p s is assumed expressed in atmospheres, h e and h w are in Btu/lb, and the 
left hand side is given in Btu/ft3/2 sec. Equation (468) is an approximation to 
the Fay- Riddell formula for hemispherical models. The other additional value 
printed is an estimate of the heat flux to a flat plate with a sharp leading edge 
at zero angle of attack: 

[q w V*] F p = — ,7 -- ( P M “l) 172 ( h r~ h w> (469) 

(N* p p 2/3 
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Here x represents the distance from the leading edge, 
Np r , p , n are evaluated at a reference temperature 



and the -gas properties 
T* defined by 

(470a) 


■h 


* 


- -L ( li e + h w ) + 0.22 (h r - h e ) 


(470b) 


h f = h e + P.r 1/2 (h s -h e ) 


(470c) 


where h* is the reference enthalpy and h r the recovery enthalpy. These formulas 
were obtained from Schlichting (ref. 21, p. 286), Jakob (ref. 69, pp. 428-430), 
and the first edition of Hayes and Probstein (ref. 25, pp. 296-297). Since they 
do not take low Reynolds number effects into account, they are expected to be- 
come inaccurate at large nozzle expansion ratios where the gas density becomes 
very low. 


8.1.5 Low-Density Effects 

The preceding stagnation point heat flux formulas are based on boundary 
layer theory and assume the existence of an inviscid shock layer. At large noz- 
zle expansion ratios, the flow becomes highly rarefied and this system of 
approximations breaks down. According to Cheng (ref. 70), the influences of 
low-density effects upon the flow in the stagnation region of a blunt body can 
be classified into the following regimes of continuum flow based upon the value 
of a parameter § K 2 : 

Regime 0: Inviscid shock layer plus boundary layer 

S> K 2 extremely large 

Regime 1 : Vorticity- interaction and viscous-layer effects 

§ K 2 > 0(1) (471) 


Regime 2: 


Merged-layer regime 


0(f p ) ^ §K 2 ^ 0(1) 


(472) 


The notation O(x) means "order of x". 
flow description becomes inappropriate. 


For § K 2 less than OU p ) , 
The parameters € and 


the continuum 
are defined by 


& m 1 (473a) 

2 y 

k 2 Pi u i R b /T* Pp \ (473b) 

H \ T 0 n* ) 

where y denotes the specific heat ratio and quantities with an asterisk super- 
script are evaluated at a reference temperature behind the shock. 
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The quantity 6 K 2 is computed by NATA to provide a basis for evaluating the 
importance of low-density effects upon the stagnation-point heat flux. Based 
on a survey of experimental data given in figure 4.1 of reference 70, the in- 
crease in q s due to vorticity-interaction and viscous- layer effects should be 
less than 10 percent for ^K 2 ^ 15, and less than 20 percent for § K 2 ,> 8. For 
S K 2 j< 1.5, the heat flux decreases with further decreases in & K 2 as a result of 
merged-layer effects. However, the heat flux should be within about 30 percent 
of the value based on boundary layer theory so long as £ K 2 ^ 0.1. 

In the calculation of §K 2 by the code, & is computed from the shock density 
ratio ip : 

S = -A_ (474) 

1 + £ P 

i* . 

which is based on the formula fpij m = (y-l)/(y + 1) for the density ratio of a strong 
shock in a perfect gas. The factor (T /xq/To #**) in (473b) is approximated by 
unity. Since the model radius Rb is not specified in the code, € K 2 is computed 
and printed for Rb = 0.3048 m (1 ft). The parameter value S K 2 for any specific 
model can be obtained by multiplying the printed value by the model radius in 
feet. 


8.1.6 Shock Standoff Distance 

NATA computes the ratio A/Rb of the shock standoff distance to the body 
radius for each type of model. This ratio is of interest because it is an ob- 
servable of model tests in wind tunnels. In general, it depends upon the shock 
strength as- measured by the density ratio e p , the model geometry, and the gas 
equation of state. The equation of state dependence is weak and is neglected in 
NATA. 


For hemispherical models, A/Rb is computed from the following curvefit to 
Van Dyke's results (ref. 71 and ref. 64, p. 462): 


A 

% 


0.78 e n 

for r < 

0.19 

P 

P - 

0.78 e p [l + 3.5 U p -0.19) 1 ' 6 ] 

f°rr p > 

0.19 


(475) 


Figure 35 compares Van Dyke's results with values obtained using the computer 
program documented in reference 63 and with Vinokur's results (refs. 65, 66) 
based on the constant-density approximation. The two sets of results from 
numerical inviscid flowfield solutions are in good agreement. Vinokur's results 
are too low by about 17 percent. 


For flat-faced cylindrical models in axisymmetric flow, A/Rb is calculated 
using a fit to Vinokur's results (refs. 65, 66), multiplied by a correction fac- 
tor of 1.17 based on the comparison for spherical models in figure 35. The 
corrected fit is 


A 

H 


1,12 


, 0.39 
P 


(476) 
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84-1806 e p = p-|/p 2 

Figure 35 SHOCK STANDOFF DISTANCE ON A SPHERE 
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For cylindrical models in two-dimensional flow, equation (6.5.14) of Hayes and 
Probstein (ref. 64) is used*: 

• 1 

1 f (477) 


Ru 




For flat-faced models in two-dimensional flow, A/R^ is calculated from the 
formula (ref. 64, p. 310) 


— = - 1-23 e* (0.79 yfr p ) 


(478) 


8.1.7 Comparisons with Experimental Data 

Scott (ref. 72) has measured pressure and heat transfer distributions over 
a 5.08-cm radius hemisphere in the NASA-Johnson Space Center 10 Mw arc tunnel. 

The test conditions are summarized in talkie I. The stagnation enthalpy was 
determined from careful energy balance and mass flow measurements. The models 
were positioned 2.5 cm downstream of the 0.508-m diameter nozzle exit. The 
0.521-m effective test section diameter listed in table I was obtained by extra- 
polating the 15-degree expansion angle of the conical nozzle from the actual 
exit plane to the model station. 

The free-stream conditions listed in the second part of table I are based 
upon runs with the current version of NATA and differ slightly from the condi- 
tions tabulated in reference 72. The initial conditions for the NATA solutions 
were the stagnation enthalpy and mass flow. Effects of the boundary layer 
displacement thickness upon the effective flow cross section at the throat were 
taken into account in the determination of the effective reservoir temperature 
and pressure from these data. The inviscid flow solution assumed chemical 
non-equilibrium. The boundary layer was included in the solution and had a 
major effect upon the effective area ratio at the test station. 

Table II summarizes the measurements of model stagnation point conditions 
in reference 72, and presents NATA results for comparison. The conditions listed 
as NATA results are based on calculations for a chemically frozen shock layer and 
a frozen boundary layer. These assumptions concerning the chemical state of the 
flow around the model are justified by the low pressure. 

The NATA predictions of stagnation pressure on the model are about 6 per- 
cent higher than the measurements in both cases. Since the stagnation pressure 
varies with the effective area ratio, this agreement is better than might have 
been expected in view of the uncertain accuracy of approximations used in the 
calculation of the displacement thickness for the boundary layer on the nozzle 
wall. 


Scott (ref. 72) measured heat transfer to models coated with both nickel 
and Teflon in order to investigate catalytic effects. The catalytic efficiency 

•Hayes and Probstein give the formula in terms of the shock radius, R s . In (477), it has been re-expressed in terms of body radius R^, 
assuming R s = R)-, + A, 
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TABLE I 


TEST CONDITIONS FOR STAGNATION CONDITION MEASUREMENTS 


Case 

1 

2 

Data 



Standard nozzle number 

10 

10 

Stagnation enthalpy, MJ/kg 

25.3 

17.6 

Total mass flow, g/sec 

13.6 

22.7 

Test section diameter, m 

0.521 

0.521 

Axial coordinate, m 

0.950 

0.950 

NATA Results 



Mach number 

10.1 

9.3 

Static temperature, °K 

267 

321 

Static pressure, atm 

3.65 x KT 5 

5.53 x 10“ 5 

Free-stream density, kg/m^ 

3.05 x 10~5 

4.41 x 10-5 

Free-stream velocity, km/sec 

4.35 

4.07 

Geometric area ratio 

356 

356 

Effective area ratio 

161 

202 

Free-stream mole fractions 



n 2 

0.270 

0.456 

N 

0.461 

0.236 

0 

0.269 

0.308 

Other species 

< 3.2 x 10 -5 

< 1.1 x 10 -5 
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TABLE II 


COMPARISON OF CALCULATED AND MEASURED 
STAGNATION CONDITIONS 


Case 

1 

2 

Measurements 



Pressure, atm 

4.83 x 1CT 3 

6.26 x 10“ 3 

Heat flux (nickel), W/cm 2 

56 

55 

Heat flux (Teflon) , W/cm 2 

36 

35 

Shock standoff, A/R^ 

0.17 

0.10 

NATA Results 



Pressure, atm 

5.17 x 10“ 3 

6.62 x 10“ 3 

Heat flux (non-catalytic) , 
W/cm 2 

35.0 

35.1 

Shock standoff, A/Rj, 

0.16 

0.14 

h D /(h e - h w ) 

0.619 

0.516 

N Le 

0.824 

0.961 

Heat flux for 0=1, W/cm 2 

85.3 

71.6 

£ K 2 

1.0 

1.3 

Rep 

130 

220 
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of Teflon for atom recombination is much lower than that of nickel. The NATA 
calculations of stagnation point heat flux to a non-catalytic surface agree with 
the data for the Teflon-coated model to within about 3 percent. The higher heat 
transfer rates for nickel could be accounted for by assuming values for $ less 
than unity in equation (466). Table II lists values of heat flux for $ = 1, 
calculated from the values for $ = 0 using (466) and the tabulated values of 

h lA h e “ h w) and N Le • 

The close agreement of the NATA stagnation point heat flux results with 
Scott's data for the Teflon-coated models has to be regarded as fortuitous. The 
low density parameter & K 2 , calculated from equations (473) and (474) , is about 
unity in both cases. Thus, the flow around the model is on the boundary of the 
merged layer regime, and low-density effects upon the heat transfer are ex- 
pected to be of substantial magnitude. 


as 


Table II includes values of Cheng's Reynolds number parameter Rep , defined 


Re t 


PI u l R b 


<p/2 


1 + 




S K 2 


(479) 


Cheng has presented a correlation of stagnation point heat transfer as a function 
of this parameter (fig. 9 in ref. 73). For Rep in the range from 100 to 200, 
the ratio q/q^L of the actual heat transfer rate to its value based on boundary 
layer theory is about 1.0 ±0.1, within the scatter of the data. At somewhat 
higher values of Rep in the vicinity of 500, q/qBL peaks at about 1.1 ±0.2, be- 
fore beginning its asymptotic approach toward unity as Rep -♦ °° . 

Scott's calculations (ref. 72), based upon a two-dimensional viscous flow- 
field model, indicate that the heat flux to the Teflon-coated model is 20 to 30 
percent higher than the flux to a completely non-catalytic surface. Thus, it 
appears likely that the NATA heat flux results are actually about 20 to 30 per- 
cent high under these low-density conditions. 


Table II also gives experimental values of the ratio A/Rp of the shock stand- 
off distance to the model radius, based upon the luminous region of the shock 
layer in photographs reproduced in reference 72. The NATA predictions of A/R b 
are roughly in agreement with the data. 


8.2 Wedge-Model Calculations 

NATA provides calculations of conditions on a blunt wedge with arbitrary 
leading edge radius and arbitrary (but not excessively large) inclination of the 
wedge surface to the flow. The principal quantities calculated are the heat 
flux, the surface pressure, the displacement thickness of the boundary layer on 
the wedge, and the shock ordinate, all as functions of distance from the leading 
edge. The positions on the wedge where these calculations are done are under 
input control. These positions can be specified in terms of an initial distance 
from the leading edge and a constant increment to generate a sequence of evenly 
spaced points, or by direct input of the positions. Both methods can be used 
together, if desired. 
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The locations in the wind tunnel at whiqh wedge-condition calculations are 
performed are controlled by the same inputs as are used to specify locations for 
stagnation-point model condition calculations. The leading edge of the wedge 
is assumed to be positioned at each of the specified locations. The actual 
variation in free stream conditions with distance along the wedge is neglected; 
in each case, the wedge is assumed to be inserted into a uniform flow with the 
velocity, density, etc. , prevailing at the location of the leading edge. 

The control variables in NATA are preset such that wedge calculations are 
not normally pefforined. To invoke the wedge calculations, it is necessary to 
specify lanfe ^r^'more 1 model' ^'s£t£pns f ^ncr "to input positive values for NANGLE ("the 
number of specified angles of attack) and NRADLE (the number of leading edge ^ ' ( 
radii). Normally, stagnation-point model condition calculations are done in 
addition to the wedge calculations. However, the stagnation-point calculations 
can be suppressed. 

The analysis, underlying the wedge condition calculations is explained ip 
the remaining parts of this section. Section 8.2.1 summarizes the classical 
theory of high-density supersonic flow over wedges. Section 8.2.2 discusses the 
low density flow over wedges including bluntness and displacement effects, as 
formulated by Cheng , et al. (ref. 74) and Kemp (ref. 75). Section 8.2.3 pre- 
sents an analytical representation of the solutions of the Cheng equation; and 
Section 8.2.4 explains certain modifications of the Cheng-Kemp results as used 
in NATA. Section 8.2.5 discusses limits to the validity of the calculations, 
and Section 8.2.6 presents a comparison of NATA wedge calculations with experi- 
mental data. 


8.2.1 Classical Wedge Theory _ * s . 

Supersonic flow over a wedge is relatively easy to calculate If the leading 
edge is sharp and the gas density is high. For a sharp leading edge, the shock 
is attached to the body at the leading edge if the Mach number is high enough. 

If the density is high, the boundary layer thickness is small in comparison with 
the shock layer thickness (except in a small region near the leading edge); thus, 
the shock is essentially straight, the density and flow velocity behind it are 
constant, and the streamlines behind the shock are parallel to the wedge sur- 
face. Under these circumstances, the conditions between the shock and the sur- 
face can be determined using oblique shock theory (ref. 64, p. 217-221 or ref. 

42, p. 85-89), The solution for the conditions behind the shock can be reduced 
to the solution of a cubic equation in the case of a perfect gas. 


8.2.2 Low-Density Hypersonic Flow Over a Wedge 

Wedge models to be tested in arc-heated wind tunnels must have blunt leading 
edges to avoid excessively high heat transfer and melting or ablation of the 
leading edge. Also, if the models are fairly large, they must be tested at 
stations where the nozzle area ratio is large, and thus where the Mach number is 
high and the free stream density low. Hypersonic low-density flow over a blunt 
wedge is dominated by effects which are neglected in the classical theory for 
sharp wedges. 
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The shock oyer a blunt leading edge is detached, and -nearly normalin the 
region near the stagnation streamline. The gas which; passes throughthis strong 
shock over the leading edge undergoes a much larger Increase In entropy than 
that which flows through the highly oblique portion of the shock further to the 
rear. This shock-heated gas forms a low-density "entropy layer" between the 
outer edge of the boundary layer and the inner boundary of a conventional thin 
shock' layer. 


At low free-stream densities, the thickness of the boundary layer on the. 
wedjge is xomparable with that of jthe inviscid shock layer. The thick boundary 
layer r.;.Spi'aces the inviscid flow awSy from the body surface, by an amount which 
increases with distances from the leading edge. This displacement results in 
an induced pressure which is highest near the leading edge, because the Shock, 
is less oblique than it would be in the absence of the boundary layer. 


Cheng et al. (ref. 74) formulated an analytical treatment of the combined 
effects of bluntness and boundary layer displacement on the hypersbnic flow over 
a wedge. Their theory treated the boundary layer using local flat-plate simi- 
larity with the simplifying assumptions of Prandtl number unity, viscosity 
proportional to absolute temperature, and constant wall temperature. The Nswton- 
Buaemann pressure law (ref. 64) was applied to the outer thin shock layer* The 
continuity equation for the entropy layer was approximated using the assumption 
y - 1«1, where y denotes the specific heat ratio for the gas. The theory was . 
found to predict the trend of experimental heat transfer data very well, but the 
predicted heat transfer was too low by about 40 percent . 

Kemp (ref. 75) modified the Cheng theory to include the effects of y 4 1 to 
first order. This modification brought the theoretical predictions into good 
agreement with experimental data, not only for heat transfer but also for sur- 
face pressure and shock shape. The wedge calculations in NATA are based on 
Kemp's analysis, which will be referred to as the Cheng-Kemp theory. 


The fundamental relation in the theory is the Cheng equation: 
d / dz\ / dz ' 1/2 




= l 


(480) 


Here 2 is a non-dimensional variable proportional to the shock ordinate, while (. 
is a non-dimensional distance from the leading edge of the wedge, and the 
constant parameter r is proportional to the angle of attack,^ The definitions 
of these quantities in terms of physical variables differ a little between the 
original and modified versions of the theory. In Kemp's form of the theory, in 
the notation of Boger and Aiello (ref. 76), 
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r = — 
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8« p (rN) 4 Y s /k 3 t 
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(48ia) 


(481b) 


(481c) 
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where 

A 

V 

"• f 

N 

N Re,t 

Y . 

k 

X 

a 

C 

T* 

T 0 

T 

w 

y 

M 



y+l 


2 



0.664 + 1.73 T^Tq 


M 



:d' ■ 


PlUjt/^i > Reynolds number based on leading-edtge thickness 


ordinate o£ shock (figure 36) 



drag coefficient of the leading edge (4/3 for a cylindrical lead- 
ing edge, 2 for a flat-faced edge) 


coordinate parallel to the free-stream velocity, zero at the 1 
leading edge (see figure 36) 


angle of attack in radians (see figure 36) 
P(T‘) . Tl_ 
m(Ti) X* 


Cheng's reference temperature, 

T* = Tq [1 + 3 (T w /Tq)]/6 

stagnation temperature 

wall temperature 

ratio of specific heats 

free-stream Mach number 

thickness of leading edge (see figure 36) 
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Figure 36 GEOMETRY OF FLOW OVER A BLUNT WEDGE 


8.2.3 Analytical Curvefit to the Solution of the Cheng Equation 

For zero angle of attack (r = 0) , equation (480) has the analytical solution* 
(ref. 74): 


z = 2 |_^ x3/ 2 - j 4 + y/\ - in (1 + 

£ = j d + A) 4 - ~ (1 + n/X) 3 +9(1 + ~ (1 + yfi) 

+ — in ( 1 + y/\) - 4>/A in (1 + yjii) + 2 [in (1 + aA)] 2 + — 

3 9 



(482a) 


(482b) 


in terms of a parameter 4 defined by 



(483) 


For T 5^ 0, the solution of (480) cannot be obtained in analytical form. However, 
the solution can be approximated by assuming that 

(zz')' = [(zz')']p = 0 + T 2 <*84) 

where the prime denotes d/d£ and [(zz')']p = o is evaluated as a function of £ 
from the solution (482), (483) for zero angle of attack: 

[(zz')'] r = 0 =i-(l+ y/J.) (485) 

For given £ , the parameter 4 must be obtained by numerical solution of equation 
(482b) and z must then be computed from (482a) . Figure 37 compares the approxi- 
mation (484), represented by the curves, with exact numerical solutions** of 
equation (480) for various values of r . The approximation (484) is most accurate 
for small r. For example, with r = 0.1 , the maximum error is about 6 percent. 

Even for large values of T , the error in (484) is less than about 10 percent ex- 
cept in the region where (zzO' is leveling off and beginning to approach its 
asymptotic value of T 2 for large £ . In this region, for larger, the exact solu- 
tion undershoots the asymptotic value and approaches it in an oscillatory manner. 
For r = 100 , the difference between (484) and the exact solution becomes as large 
as a factor of about 1.45 in this region. 

However, Cheng et al. (ref. 74) suggest that the oscillations in the exact 
solution, which are responsible for these large differences, are unphysical 
artifacts of the system of approximations upon which the theory is based. If so, 
the approximate formula (484) might have higher physical accuracy than the exact 
solution. Comparisons with experimental data, presented below, support this 
conjecture. 


*The equations in reference 74 corresponding to (482) contain minor errors. The correct forms were given by Cheng, H.K., et al (ref. 77) 

**These numerical solutions were computed and provided by R. C. Boger, Avco Systems Division. 
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The physical conditions at the surface of the wedge are given by the Cheng- 
Kemp theory in the form (ref. 76): 

Pw a 5 W M 2 (rN) 4 

• (=)' 


Pi y2 


> 


*6 


C H = 


y 1/2 f 2 r 5 N 6 (zz/) ' 


8 * = 


r 


5/2 k 3 


(486) 

(487) :. 

(488) l«4 


A 5 8€ p (rN) 4 (22 ' ) ' 

Hete p^/pi is the 'ratio of the pressure on the wedge surface to the free- 
stream pressure, Ch is the surface heat transfer coefficient. 


C H = 


Pl UjCho -h w ) 


(489) 


and f is the boundary layer displacement thickness. The shock ordinate Y s may 
be obtained from (481a): 


Y s = 


y k3 


(490) 


8 A 4 (rN) 4 


.7 The quantity ( zz ') ' appearing in (486) is given directly by the approximate 
solution (484) of the Cheng equation. The quantity zz' in (487) and (488) can; 
be obtained from (484) and (483): s f> 


/ ; 


;zz' = | t(zz')/] 

'0 


r = 0 


d£ + T 2 C 


= [zz']^ +. r 2 c= \+r 2 c 


(491) 


The variable z in (490) is then given by a further integration: 


jz 2 = j [zz' 


] d<+ — r 2 1 2 

r=o 2 


z 2 = iz 2 ] + r 2 <r 2 
r=o 


(492) 


where [z 2 ]* ^ may be obtained by squaring equation (482a). 

Figures 38 and 39 compare the approximate solution based on (484) with ex- 
perimental data on the pressure and heat transfer distributions over the surface 
of a very blunt wedge in a hypersonic air stream. The experimental points are 
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Figure 38 COMPARISON OF ANALYTICAL CURVEFIT TO THE CHENG EQUATION 
SOLUTION WITH EXPERIMENTAL PRESSURE DATA 




Figure 39 COMPARISON OF ANALYTICAL CURVEFIT TO THE CHENG EQUATION 
SOLUTION WITH EXPERIMENTAL HEAT TRANSFER DATA 


from figure 10 of Kemp's paper (ref. 75). For a flow which is dominated by the 
effects of leading-edge bluntness, the non-dimensional coordinate C is small, 
and the left-hand si de o f the Cheng equation (480) is dominated by its first 
term. If the term 'Jzz' is neglected, the solution of the Cheng equation for 
T = 0 is found to be 


z = ( 9 / 2) 1/3 £ 2/3 

Then (484) gives for general F 

{zt'r = (2/9) 1/3 <r~ 2/3 + r 2 

and 

22 ' = 6 1/3 c 1/3 + r 2 £ 


(493) 


(494) 


(495) 


If these expressions are substituted into (486) and (487) , there results after 
some algebraic manipulation 

2 


y +1 y M 2 a 2 Pi 


1 p w 0.382 

= 1 + 


y+ 1 Vy 


1. M 2 I e p k 


-H 


/ 0.382 \ 
0 ,32 j 1 *—) 

t 1/4 y/ 1.145 +£ 




(496a) 


(496b) 


•‘*1 


where 


a 2 (ill • 

\ 2 V kt / 


Xt M 3 Vc/N Re> ( 


2/3 


(497a) 

(497b) 


The unbroken curves in figures 38 and 39 represent equations (496a) and (496b) , 
respectively, The agreement with the experimental data is similar to that shown 
in Kemp's original figure 10. In the region £— 1 , where the pressure curve in 
figure 38 is levelling out, the agreement in both pressure and heat flux is bet- 
ter than is obtained with the exact solution of the Cheng equation (shown by the 
dashed curves in figures 38 and 39) . 


8.2.4 Modifications of the Cheng-Kemp Formulas 

Since (zz 0 ' is proportional to the pressure and C is proportional to the 
distance from the leading edge, figure 37 may be regarded as a non-dimensional 
plot of the pressure distribution over the surface of the wedge. The parameter T 
corresponds to the angle of attack. For F = 0 , the pressure falls as C~ 2 ^ for 
very small C and as C~^' 2 for very large £ . For F non-zero and positive, the 
pressure levels off to an asymptotic value, corresponding physically to the 
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classical wedge pressure, 
sure is given by 

= Ay M 2 a 2 

X oo 



From (484) , (486) , and (481c) the asymptotic pres- 

(498) 


For small angles of attack, classical oblique shock theory (ref. 64, p. 218-220) 
gives the pressure on the wedge as 



1 + 


y M 2 a 2 


1 -e 


P, os 


(499) 


where e pt osis the density ratio across the oblique shock wave. 


Since A = (y + l)/2 , equations (498) and (499) agree under the following conditions 


y M 2 a 2 

> > 

l 

(500a) 

1-f 




P, os 

y- 

l 

(500b) 

e p, os = 

y + 

T 



The second of these equations implies that the shock is strong even though 
oblique. For small wedge angles of a few degrees, these conditions are satis- 
fied only for considerably higher Mach numbers than are likely to be en- 
countered in the use of NATA to correlate arc tunnel data on wedges. 


Possible modifications of the Cheng-Kemp results for the conditions on 
wedge models have been studied in an effort to obtain formulas with a wider 
region of applicability. For theoretical reasons it appears inadvisable to 
change Kemp's choice of the value (y+l)/2 for the parameter A . An alternative 
and better way to obtain the correct asymptotic wedge pressure is to modify the 
definition of the parameter T. In place of equation (481c), F is now redefined 
as 



4 A 5 (p 2 M 2 (rN) 4 


(501) 


(502) 


and (p^/Ppos is t * ie surface pressure ratio based on oblique shock theory. In 
addition, the pressure formula (486) is modified to 


Pw (z zV (503) 

— = — — - + 1 
Pi Q 

which is simply (486) with a term of unity aiided to the right-hand side. In 
the asymptotic constant-pressure region, where (zz')' is equal to T 2 , (501) and 

(503) give P^/Pl = (Pw/Pl^os • Under conditions such that equations (500) are 
satisfied, (501) and (503) reduce approximately to the original Cheng-Kemp 
formulas . 
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In the Cheng-Kemp theory, the heat transfer coefficient and displacement 
thickness are calculated from formulas involving the surface pressure. In the 
notation of reference 76, 



(504a) 


(504b) 


Thus, consistency requires that the modification (503) in the pressure formula „ 
(486) be accompanied by corresponding modifications in the formulas for Cyj and S 


Substitution of (503) and (481b) into (504) gives 
2.656 A 6 t p 2 (rN) 5 N (zz ') ' + 0 
y 2 zz' + Q £ 

y 2 k^ c \J zz ' + O C 

8A 5 e p (rN) 4 (zz')' + Q 


According to the Cheng-Kemp theory, the shock ordinate is given by 

,2 


Y s = A 


y f kt 

* ' r 

x can a + S + 


2 pVpi 


(505a) 

(505b) 


(506) 


In the original form of the theory, this relation is equivalent to the Cheng 
equation (480). If Pw/pi and S are now assumed to be given by (503) and (505b) 
in place of (486) and (488), a new formula for Y s is obtained: 


y 2 k3 t 

Y s = Ax tana + 

8 A 4 e p (rN) 4 


1 + yjzz '+ O £ 
(zz ') ' + n 


(507) 


The first term on the right can no longer be expressed in terms of T, because Y 
as defined by (501) is not simply proportional to the wedge angle a . In spite 
of the modifications, (507) does not give the correct asymptotic shock angle at 
large distances from the leading edge of a sharp plate, because the coefficient A 
is related to the density ratio for an infinitely strong shock rather than that 
for the actual oblique shock, which can be rather weak. An ad hoc further 
modification could be made to force (507) to give the shock angle in the classi- 
cal wedge limit correctly, but such a change would have an adverse effect on the 
accuracy in cases in. which the bluntness and displacement effects are important. 
Accordingly, (507) is used as it stands. 
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The modified formulas (503), (505), and (507) are no longer consistent with 
the Cheng equation (480). Thus, these formulas do not represent a new form of 
the Cheng-Kemp theory, but only a set of ad hoc modifications of its results. 


8.2.5 Limits of Validity 

The original Cheng-Kemp theory is valid in the strong-interaction regime 
under conditions such that merging and slip effects are unimportant. Kemp 
(ref. 75) gives the following criteria: 


< 0.15 (508) 

N Re, x 


X e = t p r M 3 

If the condition (508) is violated, merging ot slip effects become significant. 
If (509) is violated, the strong- interaction approximation breaks down. These 
criteria define lower and upper values of the distance from the leading edge, 
bracketing the region in which the original Cheng-Kemp theory is expected to be 
valid. The lower limit xg is determined by the onset of merging or slip effects 
from (508) , 

44M 2 C 

x» _ metres (510a) 

^ Re/m 


N 


> 1 


Re, x 


(509) 



where Re/m is the free-stream Reynolds number per metre. The upper limit x u is 
the position at which the strong interaction approximation breaks down; from 
(509), „ , , 

f 2 t 2 C 


Re/n 


metres 


(510b) 


If the free stream Mach number M is too low, these limits come together and the 
original Cheng-Kemp theory is not valid anywhere on the wedge. If xg is set 
equal to x u , equations (510) give 

Mmin - ~ 7 (51D 

vy 

The original Cheng-Kemp theory has no region of validity for Mach numbers lower 
than M m i n . The modifications to the results of the Cheng-Kemp theory discussed 
in Section 8.2.4 were designed to extend the region of applicability of the 
results beyond the point at which strong interaction theory breaks down. To 
the extent that these modifications prove to be successful, the modified formu- 
las are not limited by the criterion (509) , but can be applied even for x > x u . 

No analogous approximation has been developed to extend the utility of the 
Cheng-Kemp theory into the region where merging and slip effects are important. 
However, according to Vidal and Bartz (refs. 78, 79), free-molecule flow theory 
gives an approximate upper limit to the heat transfer to the surface of a wedge. 
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A calculation of the free-molecule heat transfer coefficient is performed in 
NATA; if the resulting heat flux is lower than that predicted by the Cheng-Kemp 
theory, the free-molecule value is printed in place of the other. From Probstein 
(ref. 80), the free-molecule heat flux is 





^ e — (S sina) + g s i na + er f (s sina))J 


1 _-(S sina) 2 


(512) 


where W is the mean molecular weight. 



(513) 


and a e is the thermal accommodation coefficient. In NATA, a e Is assumed equal 
to unity. 

Unfortunately, the pressure predicted by free molecule theory is not an 
approximate upper limit to the actual pressure on the wedge (ref. 79). In some 
cases, the measured pressure is nearly an order of magnitude larger than the 
free molecule value. Thus, the latter does not provide a criterion for limiting 
the calculated pressures. 


8.2.6 Comparison with Experimental Data 

Scott (ref. 81) has measured the heat flux to a blunt wedge in the NASA 
Johnson Space Center 10 1*7 Arc Tunnel Facility. The test conditions are summa- 
rized in Table III. The stagnation enthalpy of the flow was determined from 
careful energy balance and mass flow measurements. The leading edge of the wedge 
model was positioned 2.5 cm downstream of the 0.635-m diameter exit orifice of 
the nozzle. The effective test section diameter given in Table III was calcu- 
lated by extrapolating the nozzle expansion to the position of the leading edge. 

The reservoir conditions were determined in the NATA solution from mass 
flow and stagnation enthalpy inputs. The non-equilibrium flow solution included 
the boundary layer. Table IV compares the NATA predictions of conditions on the 
wedge with Scott's experimental measurements of the heat flux to the surface of 
the wedge. Measurements of surface pressure were not performed. In Table IV, x w 
denotes the distance along the surface of the wedge from its leading edge. The 
values tabulated under "Cheng-Kemp Modified" are based on the modified Cheng- 
Kemp formulas (501) , (502) , (503) and (505) . Those given under "Cheng-Kemp Un- 
modified" are based oh the original Cheng-Kemp theory, equations (486)- (488). 

The heat flux calculated from the modified Cheng-Kemp formulas is about 20 per- 
cent low. The original Cheng-Kemp theory gives fluxes too low by 26 percent. 
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In this case, the condition (510b) for validity of the strong interaction ap- 
proximation is satisfied for x w up to 48 cm. Thus, the unmodified Cheng-Kemp 
theory is applicable over the entire region covered by the experimental data, 
so that the differences between the results from the modified and unmodified 
formulas are small. The condition (510a) predicts that merging effects become 
significant in this case for x w less than 14 cm. 
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TABLE III 


TEST CONDITIONS FOR WEDGE HEAT TRANSFER MEASUREMENTS 


Data 

Standard nozzle number 
Stagnation enthalpy, MJ/kg 
Total mass flow, g/sec 
Test section diameter, m 
Axial coordinate, m 
Pitot pressure, atm 
Wedge angle of attack, degrees 
Wedge leading edge radius, cm 
NATA Results 
Mach number 
Static pressure, atm 
Pitot pressure, atm 
Free-stream density, kg/m^ 
Free-stream velocity, km/ sec 
Free-stream temperature, °K 
Frozen specific heat ratio ( y ) 
Geometric area ratio 
Effective area ratio 
Free-stream mole fractions 

n 2 

N 

0 

Other species 


10 

13.9 

90.7 

0.648 

1.174 

1.43 x 10~ 2 
15 

0.95 

9.7 

1.09 x lO" 4 
1.40 x 10" 2 
9.26 x 10" 5 

4.07 
329 

1.471 

551 

393 

0.584 
0.081 
0.335 
< 1 x 10-4 
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TABLE IV 


CALCULATED AND MEASURED CONDITIONS ON A WEDGE MODEL 


Quantity 

Cheng-Kemp 

Modified 

Cheng-Kemp 

Unmodified 

Measured 

Angle of attack parameter, T 

• ' Til. 6 

11.2 

:!* .» T f *■ :f » 

Conditions at x w «■ 18.8 cm 




P w , atm 

2.24 x 10" 3 

2.03 x 10~ 3 

- 

q w , W/cm 2 

8.21 

7.66 

10.44 

S , cm 

1.05 

1.12 

- 

Conditions at x w = 26.4 cm 




P w , atm 

2.10 x 10 -3 

1.89 x 10 -3 

- 

% , W/cm 2 

6.88 

6.40 

8.62 

8* , cm 

1.25 

1.34 

- 

Conditions at x w = 34.1 cm 




p w > atm 

2.01 x 10 -3 

1.80 x 10 -3 

- 

q w , W/cm 2 

6.05 

5.62 

7.49 

8 , cm 

1.42 

1.53 

- 
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